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Abstract 


A  Bayesian  probabilistic  inferential  framework  provides  a  natural  and  logically  consistent 
method  for  source  reconstruction  that  fully  utilizes  the  information  provided  by  a  limited 
number  of  noisy  concentration  data  obtained  from  a  network  (or,  array)  of  detectors.  This 
report  addresses  the  application  of  this  framework  to  the  difficult  problem  of  estimating 
the  parameters  of  an  a  priori  unknown  number  of  sources,  using  an  array  of  detectors. 
To  this  purpose,  Bayesian  probability  theory  is  used  to  formulate  the  full  joint  posterior 
probability  density  function  for  the  number  of  sources  and  the  parameters  (e.g.,  location, 
emission  rate,  activation  and  deactivation  times)  that  describe  each  source.  A  simulated 
annealing  algorithm,  applied  in  conjunction  with  a  reversible-jump  Markov  chain  Monte 
Carlo  technique,  is  used  to  draw  random  samples  from  the  posterior  probability  density 
function.  By  calculating  the  marginal  posterior  probability  distribution  of  the  number  of 
sources  from  these  samples,  a  maximum  a  posteriori  estimate  Ns  for  the  number  of  sources 
can  be  obtained,  and  all  samples  of  source  distribution  models  with  exactly  Ns  discrete 
sources  can  be  used  to  provide  best  estimates  for  the  source  parameters  (along  with  their 
associated  uncertainties).  The  method  is  validated  against  a  real  dispersion  experiment 
involving  various  combinations  of  multiple  source  releases  conducted  under  a  multinational 
cooperative  FUsing  Sensor  Information  from  Observing  Networks  (FUSION)  Field  Trial 
2007  (FFT-07)  undertaken  at  US  Army  Dugway  Proving  Ground  (DPG)  in  September 
2007. 
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Resume 


Le  schema  probabiliste  inferentiel  bayesien  est  une  methode  naturelle  et  logiquement  con- 
stante  de  reconstruction  des  sources  qui  utilise  au  maximum  l’information  procuree  par 
un  nornbre  lirnite  de  donnees  de  concentration  obtenues  d’un  reseau  de  detecteurs.  Ce 
rapport  traite  d’appliquer  ce  schema  au  probleme  ardu  de  l’estimation  des  parametres  d’un 
nornbre  a  priori  inconnu  de  sources  a  l’aide  d’un  reseau  de  detecteurs.  C’est  dans  ce  but 
qu’on  utilise  la  theorie  bayesienne  des  probabilites  pour  formuler  la  fonction  de  densite 
conjointe  a  posteriori  pour  le  nornbre  de  sources  et  les  parametres  (ex.  :  endroit,  taux 
d’emission,  heures  d’activation  et  de  desactivation)  qui  decrivent  chaque  source.  On  utilise 
un  algorithme  recuit  sirnule  applique  en  conjonction  avec  une  technique  Monte  Carlo  par 
chaine  de  Markov  a  sauts  reversibles  pour  tirer  des  echantillons  au  hasard  au  rnoyen  de  la 
fonction  de  probability  a  posteriori.  En  calculant  la  distribution  de  probability  a  posteriori 
du  nornbre  de  sources  a  partir  de  ces  echantillons,  on  peut  obtenir  V estimation  du  maximum 
a  posteriori  Ns  du  nornbre  de  sources  et  on  peut  utiliser  tous  les  echantillons  provenant  des 
modeles  de  distribution  de  sources  avec  exactement  Ns  sources  discretes  pour  procurer  les 
meilleures  estimations  concernant  les  parametres  des  sources  (ainsi  que  les  fluctuations  qui  y 
sont  reliees).  La  methode  a  ete  validee  par  une  experience  consistant  en  une  dispersion  reelle 
comprenant  plusieurs  combinaisons  d’emissions  de  sources  multiples  ;  cette  experience  a  ete 
conduite  au  rnoyens  d’essais  pratiques  en  2007  (FFT-07)  par  les  reseaux  d’observation  de 
capteurs  d’information  FUSION  (FU sing  Sensor  Information  from  Observing  N etworks)  et 
effectues  sur  le  polygone  d’essais  de  l’arrne  des  Etats-Unis,  Army  Dugway  Proving  Ground 
(DPG),  en  septembre  2007. 
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Executive  summary 


Validation  of  a  Sensor-Driven  Modeling  Paradigm  for  Multiple 
Source  Reconstruction  with  FFT-07  Data 

E.  Yee;  DRDC  Suffield  TR  2009-040;  Defence  R&D  Canada  -  Suffield;  May  2009. 

Background:  A  critical  capability  gap  in  current  emergency  and  retrospective  manage¬ 
ment  efforts  directed  at  terrorist  incidents  involving  the  covert  release  of  chemical,  biological 
and  radiological  (CBR)  agents  into  the  atmosphere  is  the  ability  to  determine  the  number 
of  sources  and,  for  each  of  these  sources,  the  location,  amount  of  agent  released,  and  the 
time  of  release,  following  the  detection  of  an  event  by  a  network  of  CBR  sensors.  This  is 
the  so-called  source  reconstruction  problem  (also  referred  to  as  the  source  characterization 
or  source  inversion  problem  in  various  studies).  In  order  to  address  this  capability  gap, 
a  multinational  collaborative  program  for  the  development  of  methodologies  for  source  re¬ 
construction,  involving  the  fusion  of  (usually  noisy)  CBR  concentration  measurements  from 
remote  and  deployable  networks  of  sensors  with  model  concentration  data  obtained  from 
advanced  atmospheric  dispersion  models,  has  been  included  as  a  task  under  The  Technical 
Cooperation  Program  (TTCP),  CBR  Defence  Group  Technical  Panel  9  (TP-9)  on  Hazard 
Assessment. 

Principal  results:  In  this  report,  we  address  the  problem  of  source  reconstruction  for  the 
difficult  case  of  multiple  sources  when  even  the  number  of  sources  is  unknown  a  priori.  A 
new  methodology  is  developed  to  solve  this  problem.  This  methodology  provides  simulta¬ 
neous  estimates  for  the  number  of  sources  and  for  the  parameters  (e.g.,  release  location, 
emission  rate)  that  characterize  each  source.  The  method  is  successfully  validated  against 
a  real  dispersion  experiment  involving  various  combinations  of  multiple  source  releases. 
This  field  experiment  was  conducted  under  a  multinational  cooperative  FUsing  Sensor 
Information  from  Observing  Networks  (FUSION)  Field  Trial  2007  (FFT-07),  which  was 
designed  and  sponsored  by  TTCP  CBR  Defence  Group  TP-9  and  conducted  at  US  Army 
Dugway  Proving  Ground  (DPG)  in  September  2007. 

Significance  of  results:  The  algorithm  described  in  this  report  can  be  used  to  interpret 
agent  concentration  measurements  obtained  from  an  array  of  CBR  detectors  in  order  to 
estimate  unknown  source  characteristics  and  quantify  uncertainties  in  this  estimation.  The 
technique  can  be  used  to  estimate  the  location,  emission  rate  and  duration  of  clandestine 
CBR  agent  release(s).  Once  the  sources  have  been  characterized  in  terms  of  their  vari¬ 
ous  descriptive  parameters  (e.g.,  release  location,  emission  rate),  this  information  can  be 
subsequently  used  with  dispersion  models  to  predict  the  future  agent  transport  in  the  at¬ 
mospheric  environment  and  to  construct  toxic  corridors  with  specified  confidence  levels  in 
order  to  support  various  aspects  of  the  decision-making  process  for  emergency  managers 
and  first  responders  (civilian  and  military).  The  comprehensive  source  reconstruction  tool 
developed  here  provides  a  general  methodology  for  the  fusion  of  CBR  detector  data  with 
model  predictions  of  agent  dispersal  in  the  atmosphere.  The  application  of  this  sensor-driven 
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modeling  paradigm  will  undoubtedly  result  in  a  more  complete  situational  awareness  in  the 
operational  CBR  environment.  Indeed,  the  sensor-driven  modeling  paradigm  developed 
herein  can  be  integrated  potentially  into  operational  warning  and  reporting  (information) 
systems  that  combine  automated  data  acquisition,  analysis,  source  reconstruction,  display 
and  distribution  of  CBR  hazard  prediction  and  associated  decision-support  products.  This 
will  undoubtedly  lead  to  a  greatly  improved  Common  Operating  Picture  (COP)  in  the  CBR 
battlespace. 

Future  work:  The  next  step  is  to  develop  a  fully  operational  implementation  of  the  source 
reconstruction  capability  described  in  this  report  and  to  incorporate  this  operational  ca¬ 
pability  into  the  integrative  multiscale  urban  modeling  system  implemented  in  the  com¬ 
putational  infrastructure  at  a  government  operations  facility  (Environmental  Emergency 
Response  Section  at  Canadian  Meteorological  Centre).  This  multiscale  urban  modeling 
system  was  developed  under  a  previous  Chemical  Biological  Radiological-Nuclear  and  Ex¬ 
plosives  Research  and  Technology  Initiative  (CRTI)  Project  02-0093RD  and  provided  an 
advanced,  high-fidelity,  validated,  state-of-the-science  modeling  system  for  the  prediction 
of  urban  flows  and  the  dispersion  of  CBR  agents  released  into  these  flows.  The  success¬ 
ful  completion  of  the  future  work  proposed  here  will  give  a  government  operations  facility 
the  necessary  tools  to  provide  a  ‘whole-of-government’  (comprehensive)  single  authoritative 
source  for  expert  quality-assured  sensor-driven  CBR  hazard  predictions  and  concomitant 
decision-support  aids,  which  will  form  the  basis  for  making  decisions  for  responding  to 
and  mitigating  hazardous  release  incidents.  These  products  can  be  used  by  emergency 
managers,  planners  and  first  responders  (civil  and  military)  in  various  federal,  provincial 
and  municipal  agencies  for  informed  response  decision  making  in  both  domestic  and  inter¬ 
national  operations,  as  well  as  for  support  to  major  events  of  national  and  international 
significance  [e.g.,  Vancouver  Winter  Olympics,  Group  of  Eight  (G8)  Summit].  Furthermore, 
this  development  is  in  direct  alignment  with  Defence  R&D  Canada’s  Science  and  Technol¬ 
ogy  Strategy  for  a  Secure  Canada  and  contributes  to  the  solution  of  a  hard  problem  in  the 
defence  and  security  domain  to  build  a  reusable  major  events  security  capability. 


IV 
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Sommaire 


Validation  of  a  Sensor-Driven  Modeling  Paradigm  for  Multiple 
Source  Reconstruction  with  FFT-07  Data 

E.  Yee  ;  DRDC  Suffield  TR  2009-040  ;  R  &  D  pour  la  defense  Canada  -  Suffield  ;  mai 

2009. 

Contexte  :  II  existe  actuellement  une  lacune  grave  en  termes  de  capacites  en  gestion  des 
urgences  et  efforts  retrospectifs  contre  les  incidents  terroristes  d’emission  invisible  d’agents 
chimiques,  biologiques  et  radiologiques  (CBR)  dans  l’atmosphere.  Cette  lacune  consiste  en 
l’incapacite  a  determiner  le  nornbre  de  sources  et  pour  chacune  de  ces  sources,  le  lieu,  la 
quantite  d’agent  ends  et  l’heure  de  l’emission,  apres  qu’un  reseau  de  capteurs  CBR  ait  ef- 
fectue  la  detection  d’un  tel  evenement.  II  s’agit  d’un  probleme  de  reconstruction  des  sources 
(appele  aussi  dans  certaines  etudes  la  caracterisation  des  sources  ou  probleme  d’inversion 
des  sources).  Pour  combler  cette  lacune,  on  a  inclus  un  programme  de  collaboration  mul¬ 
tinational  a  la  tache  du  Programme  de  cooperation  technique,  Panel  9  (TP-9)  du  groupe 
technique  de  defense  CBR  devaluation  des  risques.  Ce  programme  consiste  en  la  fusion  de 
mesures  de  concentrations  CBR  (habituellement  bruyantes)  provenant  de  reseaux  de  cap¬ 
teurs  eloignes  et  pouvant  etre  deployes  dont  les  modeles  de  donnees  de  concentration  sont 
obtenus  a  partir  de  modeles  perfectionnes  de  dispersion  atmospherique. 

Resultats  principaux  :  Nous  traitons,  dans  ce  rapport,  de  ce  probleme  de  reconstruction 
des  sources  dans  les  cas  difficiles  de  sources  multiples  dont  le  nornbre  est  inconnu  a  priori. 
On  a  developpe  une  nouvelle  methodologie  pour  resoudre  ce  probleme ;  cette  methodologie 
procure  des  estimations  simultanees  du  nornbre  de  sources  et  de  ses  parametres  (ex.  :  le 
lieu  d’emission,  le  taux  d’emission)  qui  caracterisent  chaque  source.  On  a  reussi  a  valider 
cette  methode  en  experimental  sur  une  dispersion  reelle  ayant  des  combinaisons  variees  de 
sources  multiples  d’emission.  Cet  experience  sur  le  terrain  a  ete  conduite  au  rnoyen  d’essais 
pratiques,  en  2007,  (FFT-07)  par  les  reseaux  d’observation  de  capteurs  d’information  FU¬ 
SION  (FU sing  Sensor  Information  from  Observing  Networks).  Ces  essais  ont  ete  conqus 
et  parraines  par  le  Groupe  TP-9  de  la  defense  CBR  du  Programme  de  cooperation  tech¬ 
nique.  Ils  ont  ete  conduits  sur  le  polygone  d’essais  de  l’arme  des  Etats-Unis,  Army  Dugway 
Proving  Ground  (DPG),  en  septembre  2007. 

Portee  des  resultats  :  On  peut  utiliser  l’algorithme  decrit  dans  ce  rapport  pour  in¬ 
terpreter  les  mesures  de  concentration  des  agents  obtenus  d’un  reseau  de  detecteurs  CBR 
en  vue  d’estimer  les  caracteristiques  inconnues  de  la  source  et  quantifier  les  fluctuations  de 
cette  estimation.  On  peut  utiliser  cette  technique  pour  estimer  le  lieu,  le  taux  et  la  duree 
d’emission  d’un  agent  CBR  clandestin  ayant  ete  ernis.  Une  fois  les  sources  caracterisees  en 
termes  de  parametres  descriptifs  varies  (ex.  :  lieux  et  taux  d’emission),  l’information  peut 
etre  subsequemment  utilisee  au  rnoyen  de  modeles  de  dispersion  visant  a  predire  le  transport 
futur  d’agents  dans  le  milieu  atmospherique  et  visant  a  construire  des  corridors  toxiques 
ayant  des  coefficients  de  confiance  en  mesure  de  soutenir  certains  aspects  du  processus  de 
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prise  de  decisions  des  gestionnaires  des  mesures  d’urgences  et  premiers  intervenants  (civils 
et  militaires).  L’outil  de  reconstruction  complete  des  sources,  mis  au  point  ici,  procure  une 
methodologie  generate  pour  la  fusion  des  donnees  de  detecteurs  CBR  au  rnoyen  de  modeles 
de  prediction  de  la  dispersion  des  agents  dans  P atmosphere.  L’application  de  ce  paradigme 
de  rnodele  a  base  de  capteurs  resultera  sans  aucun  doute  en  une  reconnaissance  situation- 
nelle  plus  complete  lors  d’operations  dans  un  environnement  CBR.  En  realite,  ce  paradigme 
de  rnodele  a  base  de  capteurs,  mis  au  point  ici,  a  le  potentiel  d’etre  integre  dans  des  systemes 
operationnels  d’avertissement  et  de  transmission  (de  l’information)  qui  combinent  l’acqui- 
sition  automatisee  des  donnees  avec  les  analyses,  la  reconstruction  des  sources,  Paffichage 
et  la  distribution  des  predictions  des  risques  CBR  et  des  produits  d’aide  a  la  decision.  Ceci 
aboutira  sans  aucun  doute  a  une  grande  amelioration  de  l’image  commune  de  la  situation 
operationnelle  (ICSO)  dans  l’espace  de  combat  CBR. 

Perspectives  d’avenir  :  La  prochaine  etape  est  de  developper  l’implementation  totale 
de  la  capacite  de  reconstruction  des  sources,  decrite  dans  ce  rapport  et  d’incorporer  cette 
capacite  operationnelle  au  systeme  integratif  de  modelisation  urbaine  a  echelles  multiples, 
implements  dans  l’infrastructure  computationnelle  des  installations  operationnelles  gouver- 
nementales  (Section  de  reponse  aux  urgences  environnementales  du  Centre  meteorologique 
canadien).  Ce  systeme  de  modelisation  urbain  a  echelles  multiples  a  deja  ete  mis  au  point 
par  le  Projet  02-0093RD  de  l’lnitiative  de  recherche  et  de  technologie  (CRTI)  chimique, 
biologique,  radiologique,  nucleaire  et  explosive ;  il  s’agit  d’un  systeme  de  modelisation 
perfectionne,  de  haute  fidelite,  valide  et  d’une  science  d’avant  garde,  visant  a  predire 
les  flots  urbains  et  la  dispersion  des  agents  CBR  ends  dans  ces  flots.  La  realisation  des 
travaux  futurs  proposes  ici,  procurera  aux  installations  operationnelles  gouvernementales, 
les  outils  necessaires  facilitant  une  source  gouvernementale  autorisee,  unique  est  globale 
(comprehensive)  qui  fera  des  predictions  d’une  qualite  assuree  des  risques  CBR  provenant 
des  capteurs.  Ces  outils  procureront  aussi  les  aides  concomitantes  qui  seront  a  la  base  des 
prises  de  decisions  ayant  trait  a  la  reponse  et  a  l’attenuation  des  incidents  d’emission.  Ces 
produits  pourront  etre  utilises  par  les  gestionnaires  des  mesures  d’urgences,  les  planifica- 
teurs  et  les  premiers  intervenants  (civils  et  militaires)  des  agences  federales,  provinciates 
et  municipales  qui  seront  alors  en  rnesure  de  prendre  des  decisions  eclairees  en  matiere 
d’operations  domestiques  et  internationales  ainsi  qu’en  soutien  a  des  evenements  irnpor- 
tants  de  portee  nationale  et  internationale  (ex.  :  les  jeux  olympiques  d’hiver  de  Vancouver, 
le  sornmet  du  Groupe  des  huit  (G8)).  De  plus,  ce  developpement  s’aligne  directement  avec  la 
Strategie  scientifique  et  technologique  pour  un  Canada  securitaire  de  R&D  pour  la  defense 
Canada  et  contribue  a  resoudre  le  probleme  difficile  en  matiere  de  defense  et  de  securite 
consistant  a  construire  une  capacite  en  matiere  de  securite  des  evenements  importants  qui 
soit  reutilisable. 
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1  Introduction 


The  development  of  increasingly  more  sophisticated  sensing  technologies  for  the  monitoring 
of  the  concentration  of  hazardous  contaminants  [e.g.,  chemical,  biological,  or  radiological 
(CBR)  agents]  released  into  the  turbulent  atmosphere  has  generated  interest  in  utilizing  this 
information  for  the  reconstruction  of  the  contaminant  source  responsible  for  the  observed 
concentration  pattern.  More  specifically,  in  public  security  applications  for  countering  ter¬ 
rorist  incidents  involving  the  covert  release  of  a  CBR  agent  in  a  densely  populated  urban 
centre,  a  critical  requirement  is  the  characterization  of  the  unknown  source(s)  following 
event  detection  by  a  network  (array)  of  CBR  sensors.  These  sensors  are  placed  at  differ¬ 
ent  points  in  space  within  a  designated  region  in  order  to  function  as  detectors/monitors 
to  provide  quantitative  measurements  of  the  concentration  of  various  air  admixtures  of 
contaminants. 

For  example,  the  Department  of  Homeland  Security  (DHS)  has  deployed  (albeit  sparse)  ar¬ 
rays  of  biological  agent  sensors  in  31  (with  plans  to  expand  to  120)  cities  across  the  United 
States  as  part  of  the  BioWatch  program  [1]  in  order  to  provide  detection  and  warning  of  a 
covert  bioterrorisnr  event.  In  the  context  of  homeland  security,  the  BioWatch  program  has 
provided  the  impetus  for  recent  research  efforts  directed  towards  the  source  reconstruction 
problem  for  determination  of  the  location,  emission  rate,  and  other  characteristics  of  un¬ 
known  source(s)  of  contamination.  In  consequence,  a  multinational  collaborative  program 
for  the  development  of  methodologies  for  source  reconstruction  has  been  included  as  a  task 
under  The  Technical  Cooperation  Program  (TTCP),  CBR  Group  Technical  Panel  9  on 
Hazard  Assessment. 


Mathematically,  the  source  reconstruction  problem  is  an  inverse  problem.  Let  c(x,  t )  denote 
the  instantaneous  concentration  at  location  x  and  time  t.  The  ensemble-mean  concentration 
C(x,  t)  =  (c(x,  t))  ((  •  )  denotes  an  ensemble-averaging  operation)  is  related  to  the  source 
density  function  S(x,t)  through  a  Volterra  integral  equation  of  the  first  kind  (exact  source- 
receptor  relation): 

C(x,t)  =  f  f  p(x,  t')S(x',  t')  dtx!  dt' ,  (1) 

J to  J R3 

where  for  simplicity  it  is  assumed  implicitly  that  at  the  (arbitrary)  initial  time  to,  C(x,  to)  = 
0.  In  Eq.  (1),  p(x,t\x',t')  =  (RK(x,  t\x’ ,  t'))  is  the  kernel  function  with  RK  determined  as 
the  fundamental  solution  of  an  advection-diffusion  equation: 


— f  +  u  •  Vi4  =  kV2Rk,  (2) 

at 

with  initial  condition  RK(x,  to|x/;  to)  =  5(x  —  x')  (where  <5(x)  is  the  Dirac  delta  function).  In 
Eq.  (2),  u  is  the  instantaneous  velocity  field  which  is  assumed  to  be  turbulent  (and,  hence, 
random  and  unpredictable)  and  k  is  the  molecular  diffusivity  which  is  taken  to  be  constant. 
The  source  reconstruction  problem  is  as  follows:  given  C(x,t)  and  a  specific  p(x,t\x'  ,t'), 
what  is  the  source  density  function  S(x',t')7 
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The  right-hand-side  (RHS)  of  Eq.  (1)  defines  a  formal  operator  G  as  follows: 

(GS)(x,  t)  =  I  I  p(x,t\x  ,t')S(x.' ,t')  dx  dt' ,  (3) 

Jt0  J R3 

which  allows  us  to  write  Eq.  (1)  in  the  symbolic  form 

GS  =  C.  (4) 

Formally,  the  source  reconstruction  problem  can  be  solved  by  constructing  the  inverse 
operator  G~l .  Unfortunately,  this  direct  mathematical  inversion  is  not  possible  for  a  number 
of  reasons.  Firstly,  there  is  no  comprehensive  theory  of  turbulence  so  an  exact  form  of 
p(x,  t\x.r ,  t')  (and  of  the  operator  G)  is  not  available.  Secondly,  as  a  pure  integral  equation  the 
problem  may  have  no  solutions  (singular  operator  G),  and  when  it  does  the  solutions  may 
not  be  uniquely  determined.  Thirdly,  C(x,f)  cannot  be  specified  completely  and  precisely 
—  measurements  of  C  are  available  only  at  a  finite  number  of  space-time  points  and  these 
measurements  are  usually  noisy  suggesting  that  even  if  G~l  was  known  exactly,  the  presence 
of  noise  in  the  concentration  data  may  introduce  instabilities  into  the  solution.  For  these 
reasons,  the  source  reconstruction  problem  using  incomplete  and  noisy  concentration  data 
is  an  ill-posed  problem. 

A  mathematical  tool  that  has  been  developed  for  treating  instabilities  that  occur  in  inverse 
problems  is  regularization  [2].  In  this  approach,  the  class  of  admissible  solutions  for  S 
is  restricted  by  imposing  further  constraints  on  the  solution  (although  frequently  these 
constraints  are  somewhat  arbitrary  and  ad-hoc).  More  specifically,  regularization  of  the 
source  reconstruction  problem  involves  finding  the  “best”  result  from  among  all  those  source 
distributions  that  agree  with  the  incomplete  and  noisy  concentration  data,  by  minimizing 
a  cost  functional  of  the  form: 


J(S)  =  \\C-GS\\2  +  \<S>(S),  (5) 

where  the  first  term  on  the  RHS  of  Eq.  (5)  is  a  measure  of  the  misfit  between  the  measured 
(noisy)  concentration  C  and  the  model  concentration  GS  (where  ||  •  ||  is  some  appro¬ 
priately  defined  norm),  $(£)  is  the  regularization  functional  that  is  used  to  impose  some 
constraint  on  the  solution,  and  A  is  a  regularizing  parameter  that  imposes  a  relative  weight 
(or,  importance)  between  the  data  and  the  constraint. 

Robertson  and  Persson  [3]  and  Robertson  and  Langner  [4]  minimized  Eq.  (5)  with  4>(5)  =  0, 
using  a  four-dimensional  (space  and  time)  variational  data  assimilation  method  with  an 
adjoint  transport  equation,  to  recover  emission  rate  profiles  for  an  experiment  with  synthetic 
data  and  for  the  European  Tracer  Experiment  (ETEX)  [5],  respectively.  Seibert  and  Stohl 
[6]  and  Seibert  [7]  use  a  regularization  functional  of  the  form  4>(S')  =  |  |Sj  |2  (which  imposes  an 
upper  bound  on  the  “energy”  of  the  source  distribution  S)  to  reconstruct  the  distribution  of 
emission  rates  for  application  to  the  radionuclide  monitoring  system  implemented  under  the 
Comprehensive  Test  Ban  Treaty  (CTBT)  [8].  Thomson  et  al  [9]  investigated  the  use  of  three 
different  regularization  functionals  ^(S)  for  source  reconstruction,  with  the  minimum  of  the 
cost  function  J(S)  obtained  using  a  random  search  algorithm  with  simulated  annealing. 
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The  three  regularization  functionals  were  an  entropy  functional,  a  smoothing  functional  in 
which  the  source  strength  in  a  grid  cell  is  compared  to  the  average  of  the  source  strengths 
in  the  neighboring  grid  cells,  and  a  functional  that  imposes  an  upper  bound  on  the  total 
emission  flux.  Bocquet  [10]  used  entropy  as  a  regularizer  (viz.,  applied  the  principle  of 
maximum  entropy  on  the  mean)  and  solved  the  resulting  constrained  minimization  of  the 
cost  function  by  using  duality  theory  to  convert  the  problem  to  a  simpler  unconstrained 
minimization  problem  for  the  Lagrangian  multipliers  (which  were  introduced  to  impose 
constraints  on  the  mean  to  ensure  that  the  model  concentration  is  consistent  with  the 
measured  concentration).  Finally,  Allen  et  al  [11]  minimized  a  misfit  measure  defined 
as  the  squared  differenced  between  the  logarithmic  measured  and  modeled  concentration 
with  <h(S')  =  0.  In  effect,  these  investigators  assumed  that  it  was  known  a  priori  that 
the  unknown  source  is  a  continuous  point  source  with  unknown  source  location  (x,  y )  and 
source  strength  Q,  and  used  a  genetic  algorithm  to  select  the  source  location  and  strength 
so  as  to  minimize  the  misfit  measure. 

The  regularization  approach  selects  a  particular  source  distribution  S  by  minimizing,  max¬ 
imizing,  or  optimizing  some  form  of  cost  functional.  Unfortunately,  the  reliability  of  the 
inferred  S  cannot  be  obtained  from  any  single  selection,  no  matter  what  optimal  properties 
in  terms  of  quality  and  utility  this  single  selection  of  S  purportedly  embodies.  To  deal 
with  uncertainty  (which  reflects  also  the  non-uniqueness  of  solutions  for  S  using  incomplete 
and  noisy  concentration  data),  it  is  necessary  to  apply  Bayesian  probability  theory  to  the 
problem,  rather  than  simply  apply  regularization.  To  this  purpose,  a  probabilistic  approach 
using  a  Bayesian  inferential  scheme  that  allowed  the  uncertainty  in  the  inference  for  S  to 
be  determined  was  developed  by  Yee  [12]  and  demonstrated  using  Project  Prairie  Grass 
data  for  dispersion  over  open  terrain.  This  methodology  was  further  developed,  refined, 
and  generalized  in  subsequent  work:  (1)  application  of  the  methodology  to  complex  envi¬ 
ronments  (dispersion  in  built-up  environments)  by  Yee  [13],  Keats  et  al  [14]  and  Chow  et  al 
[15];  (2)  generalization  of  the  methodology  to  deal  with  a  non-conservative  scalar  by  Keats 
et  al  [16];  and,  (3)  application  of  the  methodology  to  source  reconstruction  for  long-range 
dispersion  on  continental  scales  by  Yee  et  al  [17].  Yee  [18]  generalized  the  methodology 
to  the  reconstruction  of  multiple  sources  when  the  number  of  sources  was  known  a  priori. 
Finally,  Yee  [19],  Yee  [20]  and  Yee  [21]  developed  the  theory  underlying  the  application  of 
a  Bayesian  probabilistic  inferential  framework  for  addressing  the  problem  of  source  recon¬ 
struction  for  the  difficult  case  of  multiple  sources  when  the  number  of  sources  is  unknown 
a  priori.  For  a  general  review  of  source  estimation  methods  for  atmospheric  dispersion,  the 
reader  is  referred  to  Rao  [22]. 

The  aim  of  this  paper  is  to  use  some  new  concentration  data,  measured  by  a  sensor  array 
consisting  of  100  detectors  for  releases  involving  multiple  sources,  to  test  the  procedure 
proposed  by  Yee  [21]  for  multiple  source  reconstruction  for  the  case  when  the  number  of 
sources  is  unknown  a  priori. 
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2  Bayesian  inference  for  source  reconstruction 


In  this  report,  we  focus  on  a  source  distribution  S  associated  with  Ns  transient  point 
sources  with  the  &:-th  source  located  at  a  vector  position  xS;fc  and  with  source  activation 
and  deactivation  times  Tb  and  Tk,  respectively,  between  which  the  source  is  releasing 
contaminant  at  a  constant  emission  rate  Qk  (k  =  1,2,...,  Ns ).  This  source  distribution  has 
the  following  explicit  form: 

N3 

S(x,  t)  =  J2  Qk*(x  -  xs,fc)  [n(t  -  Tbk )  -  H(t  -  Tek)\ ,  (6) 

fe=i 


where  TL(  ■  )  is  the  Heaviside  unit  step  function.  Now,  let  us  assemble  the  parameters 
for  this  particular  source  distribution  (consisting  of  Ns  transient  point  sources)  into  the 
following  source  parameter  vector: 


0Ns  =  (^llT^T^Q1,...,^Na,TbNs,TeN3,QNs)  G 


p6iVs 


(7) 


Furthermore,  let  0  =  ( Ns,8ns )• 


The  measured  mean  concentration  ■(»)  obtained  by  a  detector  at  receptor  location  x<^ 

and  time  is  assumed  be  the  sum  of  a  modeled  mean  concentration  signal  C (x^ ,  ;  0) 

and  noise 

di J«)  =  C^Xdj ,  ;  ©)  +  i  =  1 , 2, . . . ,  Nd  and  j  =  1 ,  2, . . .  ,  JVtW ,  (8) 


(j) 

where  1V^  is  the  number  of  detectors,  Nb  '  is  the  number  of  concentration  time  samples 
measured  at  the  i-th  detector  and 


C'(xdi>4*-’0)  =  [  dt  [  dx  C(x,t)/i(x,t|xd.,4*))  =  «C’|/i))(xd.,tJh, 

J  Jto  J©cR3  J  J 


(9) 


is  the  expected  (mean)  concentration  “seen”  by  the  detector  at  location  Xrf.;  and  time  t^, 
V  x  [to,T]  is  the  space-time  volume  enclosing  the  source  distribution  S  and  the  detectors, 
and  C*(x,  t)  is  the  ensemble-mean  concentration  determined  in  accordance  to  Eq.  (1)  for 
the  source  distribution  S  given  by  Eq.  (6).  In  Eq.  (9),  /r(x,  t|x^,  td)  is  the  spatial-temporal 
filtering  function  (of  x  and  t )  for  the  mean  concentration  measurement  made  by  the  detector 
at  location  x^  and  time  td  with 


dx  /i(x,t|xd,td)  =  1. 


(10) 


Note  from  Eq.  (9)  that  C  can  be  determined  as  the  inner  (or  scalar)  product  ((C\h))  of  the 
mean  concentration  C  and  the  detector  spatial-temporal  filtering  function  h. 


To  simplify  the  notation,  we  rewrite  the  measurement  model  of  Eq.  (8)  as  follows: 


dj  =  Cj(Q)  +  ej,  J  =  1,  2, . . .  ,  N, 


(11) 
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where  N  =  is  the  total  number  of  measured  concentration  data  and  Cj(0)  = 

C(x^.,t^;0).  In  Eq.  (11),  the  index  J  is  used  to  denote  the  label  (ordered  in  some 

regular  or  convenient  manner).  With  this  background,  the  problem  of  source  reconstruction 
reduces  to  the  following  problem:  estimate  Ns  (number  of  sources)  and  0ns  (parameters  for 
each  source)  or,  equivalently,  estimate  0  given  the  concentration  data  D  =  (di,  d?, . . . ,  d/v). 


Yee  [21]  developed  a  new  method  based  on  Bayesian  probability  theory  for  the  inference 
of  0.  The  components  of  the  Bayesian  inference  scheme  for  source  reconstruction  are 
shown  in  Figure  1.  Bayesian  probability  theory  can  be  derived  from  more  fundamental 
principles  starting  with  the  formulation  of  a  small  number  of  requirements  that  any  theory 
of  plausibility  (or  inference)  ought  to  verify.  These  requirements  (desiderata)  were  first 
provided  by  Cox  [23] ,  with  an  eloquent  description  of  the  complete  development  described 
by  Jaynes  [24]  in  his  definitive  treatise.  Within  the  context  of  the  source  reconstruction 
problem  formulated  above,  Bayes’  theorem  yields  the  following  result: 


P(0|D,J) 


p(0|J)p(Dl0,J) 
P(  D|J) 


(12) 


where  7  is  the  background  (contextual)  information  available  in  the  problem  (e.g.,  model 
that  defines  the  mapping  from  a  source  distribution  S  to  the  concentration  C .  background 
meteorology).  The  various  factors  that  appear  in  Eq.  (12)  have  the  following  interpretation. 
Firstly,  p(0|7)  is  the  prior  probability  density  function  (PDF)  for  a  proposition  0  about 
the  source,  predicated  on  the  contextual  information  specified  by  7,  with  “|”  denoting 
“conditional  upon”.  The  prior  PDF  encodes  all  the  prior  information  about  the  source 
before  receipt  of  the  concentration  data  D.  Secondly,  p(D|0,7)  is  the  likelihood  function 
and  is  the  probability  that  we  observe  the  concentration  data  D,  when  0  is  known  exactly 
(viz.,  the  source  distribution  is  known).  Thirdly,  p(D|7)  is  referred  to  as  the  evidence  and, 
in  our  case  here,  is  simply  a  normalization  constant.  Finally,  p(0|D,7)  is  the  posterior 
PDF  for  the  proposition  0  about  the  source,  in  light  of  the  new  information  introduced 
through  the  newly  acquired  concentration  data  D. 


We  are  seeking  the  posterior  PDF  p(0|D,  7),  which  encodes  our  inferences  about  0.  Because 
p(D|7)  is  simply  a  normalization  constant,  the  posterior  PDF  of  interest  in  Eq.  (12)  can  be 
specified  within  a  normalization  constant  as 

p(0|D,7)ocp(0|7)p(D|0,7).  (13) 

The  problem  now  reduces  to  the  assignment  of  p(0|7)  (prior  distribution)  and  p(D|0,7) 
(likelihood  function) . 


For  the  prior  distribution,  the  logical  independence  of  the  source  parameters  is  assumed,  in 
which  case  p(0|7)  factorizes  as  follows: 

Ns 

p(0\I)=p{Ns,9Ns\l)  =p(Ns\I)  p(Qk\I)p{xs^\I)p(T^\I)p(Tg\T^ , 7).  (14) 

k= 1 

The  prior  on  the  number  of  sources  p(Ns\I )  is  chosen  to  be  a  binomial  distribution  with 
parameter  p*  £  [0,1]  (binomial  rate),  and  with  a  domain  of  definition  between  NS)U1[n 
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p(0ID,/)oc  p(0|/)p(DI0,7) 
inference  oc  prior  x  likelihood 


Figure  1:  The  components  of  the  Bayesian  inference  scheme  for  source  reconstruction. 


(minimum  number  of  sources)  and  NSjiaax  (maximum  number  of  sources);  viz.,  p(Ns\I)  ~ 
B(p*\  Nst.m\n,NXma.Tc)1  where  the  symbol  denotes  “is  distributed  as”.  The  prior  on 
emission  rate  p(Qk\I )  for  k  =  1,2 . . .  ,NS  is  chosen  to  be  a  Bernoulli- uniform  mixture: 
p(Qk\I)  ~  BU('y,Qmax)  where  7  is  the  probability  that  the  source  is  turned  on  (viz., 
Pr{Qfc  >  0}  =  7)  and  Qmax  is  an  a  priori  upper  bound  on  the  expected  emission  rate.2 
The  prior  on  the  source  location  p(xStk\I)  for  k  =  1,2 . . .  ,  Ns  is  chosen  to  be  uniform 
(flat)  over  the  some  spatial  region  T>  C  R3,  so  p(xS)fc|7)  ~  U(T>).  Finally,  the  priors 
on  the  source  activation  (on)  and  deactivation  (off)  times  for  the  k- th  source  are  chosen 
to  be  uniform  over  [to,  Tmax]  and  [Tfefc,Tmax],  respectively,  where  Tmax  is  an  upper  bound 
on  the  time  at  which  the  source  was  turned  on  or  off  (viz.,  p(T^\I)  ~  U([to,Tmax\)  and 
p(Tefc|Tbfc ,/)  ~  U(lTfr\  Tmax])  for  k  =  1,2,...  ,NS).  Note  that  the  domain  of  definition  for 


lrThe  binomial  distribution  B(p*;  NStU1  in,  Ais,max)  has  a  probability  distribution  function  defined  as  follows: 


p(Ns\I) 


(Ns 

,max  -^Vs,  min)- 

(iVs-N3„min)!(7Vs,max-7Vs)!: 


’(l-PT 


for  N3  =  Ns^in,  NB:lnin  +  1 , . . . ,  N/s.max -  Note  that  in  this  definition,  the  standard  form  of  the  binomial 
distribution  has  been  offset  by  the  minimum  number  of  sources  AiSimin. 

2More  specifically,  a  Bernoulli-uniform  mixture  model  for  Qk  has  the  following  form: 


p(Qk\I)  =  (1  -  l)S{Qk)  +7I(o,Q*)(Qfc)/Q*, 

(A;  =  1,2,...,  Na)  and  IaI*)  is  the  indicator  function  for  set  A,  with  Ia(®)  =  1  if  x  G  A  and  Ia(*)  =  0  if 
x  A. 


6 


DRDC  Suffield  TR  2009-040 


the  uniform  distribution  assigned  to  p(T^\T^ ,  I)  explicitly  encodes  the  fact  that  the  time 
the  A;- th  source  is  turned  off  must  occur  after  it  has  been  turned  on. 


The  likelihood  function  p(D|0, 1)  encodes  all  the  information  provided  by  the  concentration 
data  about  the  unknown  source  distribution  S.  The  “noise”  component  ej  in  Eq.  (11) 
represents  the  difference  (residual)  between  the  measured  and  modeled  (or,  predicted)  mean 
concentration.  It  is  assumed  that  the  only  information  we  have  about  the  noise  component 
ej  is  that  it  has  a  known  noise  power  (or  variance)  crj.  The  variance  crj  of  this  noise  can 
be  decomposed  into  four  basic  contributions  as  discussed  by  Rao  [25]  (see  Figure  1):  (1) 
model  errors  arising  from  uncertainties  in  the  representation  of  various  physical  processes 
(e.g.,  turbulent  diffusion)  in  the  dispersion  model  used  to  predict  C*(x,  f);  (2)  input  error 
arising  from  uncertainties  in  the  values  of  model  parameters  and/or  in  the  specification  of 
the  meteorology  used  to  “drive”  the  dispersion  model;  (3)  stochastic  uncertainty  arising 
from  the  turbulent  nature  of  the  atmosphere;  and,  (4)  measurement  noise  inherent  in  the 
concentration  detector.  In  spite  of  the  complexity  of  the  noise  structure,  application  of  the 
principle  of  maximum  entropy  (see  Jaynes  [24])  to  our  state  of  knowledge  concerning  the 
noise,  results  in  the  following  Gaussian  form  for  the  likelihood  function: 


where 


P(D|0,J) 


n?=1  yfaoj  exp 


dj  -  Cj(G) 
oj 


2 


(15) 


(16) 


Using  the  forms  for  the  prior  distribution  and  the  likelihood  function  assigned  above,  the 
posterior  distribution  of  Eq.  (13)  can  be  written  as  follows: 

P(0|D,  I)  = 

oc 


The  posterior  PDF  p(@|D,  I)  embodies  the  state  of  knowledge  about  the  source  parameters 
given  the  prior  information  encoded  in  p(Q\I)  and  the  newly  acquired  concentration  data 
D,  the  latter  of  which  modulates  our  prior  belief  about  0  through  the  likelihood  function 
p(D|0,/).  In  consequence,  p(0|D,7)  allows  us  to  estimate  all  the  interesting  statistics 
about  the  source  parameter  0.  More  specifically,  the  posterior  distribution  p(0|D,  I)  allows 


p(Ns,0ns  |D,  I) 

1 


nil  v/2vrcrJ 


exp 


1  ( dj  —  Cj(@)\ 


E 

j= i 


Vj 


J 


w  (Ns, max  Akinin)!  *(jVs-jV3,min)  n  n*\Ns,max-Ns 

(N  -N  ■  V(N  -NVP  '  P  ’ 

Vivs  iVs„min/‘ViVs,max  iy,sJ- 


Ns 


Xll  C1  -l)S(Qk)  +^(0,Qmax){Qk)/Qr 


k=  1 


I trpk  rp  )(Tg) 

xiP(X,,fc)i(t0,Tmax)mfc)  1  . 

1 '  max  -^b  ) 


(17) 
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one  to  summarize  its  information  through  a  few  statistics  such  as  location  and  dispersion 
measures  which  provide  quantitative  information  on  the  central  (or,  “best”)  value  of  a 
source  parameter  and  the  uncertainty  associated  with  this  value. 

For  example,  the  maximum  a  posteriori  estimate  for  the  number  of  (discrete)  sources  can 
be  obtained  from 

Ns  =  argnrax  p(Ns\D,I),  (18) 

Na 

where  p(7Vs|D,  /)  is  the  (marginal)  posterior  probability  for  the  number  of  sources,  condi¬ 
tioned  on  the  concentration  data  and  the  background  information.  Given  an  estimate  Ns 
for  the  number  of  sources,  posterior  quantities  of  interest  about  9 ^  involve  expectation 
values  of  the  following  form: 

Wvs)>  =  f  Wfi')p{Na,eR'\DtI)dBfi'.  (19) 

In  particular,  the  choice  )  =  Qft  gives  the  posterior  mean  of  the  source  parameters, 

which  can  be  used  as  an  estimate  for  these  parameters;  viz.,  9 ^  =  ( 9 ^  ). 

Similarly,  the  choice  )  =  (8^  —  ))2,  gives  the  posterior  variance  of  the  source 

parameters,  whose  square  root  (or,  posterior  standard  deviation)  can  be  used  as  a  measure 
of  uncertainty  in  the  estimate  for  these  parameters;  viz.,  o{9^  )  =  ((9^  —  (9^ 
Alternatively,  a  p%  highest  posterior  distribution  (HPD)  interval  that  encloses  a  source 
parameter  with  p%  probability,  and  constructed  so  that  the  lower  and  upper  bounds  of 
the  specified  interval  are  such  that  the  probability  density  function  within  the  interval  is 
everywhere  larger  than  outside  it,  can  be  used  as  a  uncertainty  specification  for  a  source 
parameter. 

3  Computational  framework 


This  section  describes  the  computational  procedures  which  are  used  for  extracting  the  source 
parameter  estimates  required  for  event  reconstruction.  In  this  report,  the  background  treat¬ 
ment  on  the  computational  methodology  is  necessarily  brief.  The  reader  is  referred  to  Yee 
et  al  [17]  and  Yee  [21]  for  a  more  complete  description  of  this  topic.  There  are  two  major 
issues  in  the  computational  framework  applied  to  Bayesian  inference  for  source  reconstruc¬ 
tion  that  need  to  be  addressed;  namely,  (1)  a  computationally  efficient  methodology  for 
the  computation  of  the  source-receptor  relationship  required  in  the  determination  of  the 
likelihood  function,  and  (2)  a  methodology  for  sampling  from  the  posterior  distribution  for 
the  source  parameters. 

3.1  Fast  computation  of  source-receptor  relationship 

To  apply  the  Bayesian  inference  methodology  to  source  reconstruction,  we  need  to  relate  the 
hypotheses  of  interest  about  the  unknown  source  distribution  (encoded  in  0)  to  the  available 
concentration  data  dj  measured  by  the  network  (or  array)  of  detectors.  This  requires  the 
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Figure  2:  Commutative  diagram  illustrating  two  mathematically  equivalent  representations 
for  the  source-receptor  relationship. 


computation  of  a  modeled  (or  predicted)  concentration  Cj  as  prescribed  by  Eqs  (9)  and 
(1).  This  is  the  source-receptor  relationship  (realized  in  terms  of  an  atmospheric  dispersion 
model)  that  describes  the  mapping  Af  sr  from  the  hypothesis  space  7i  of  source  distributions 
to  the  concentration  sample  space  SN  of  concentration  data,  so  :  0  G  H  ->  D  £  SN 

(N  is  the  number  of  concentration  data). 

The  likelihood  function  given  in  Eqs  (15)  and  (16)  is  not  a  closed-form  expression  and  its 
evaluation  is  computationally  expensive  owing  to  the  fact  that  Cj  (J  =  1,2, .. .  ,N)  needs 
to  be  determined  for  a  given  source  distribution  0  [using  Eqs  (9)  and  (1)].  Moreover,  a 
simulation-based  posterior  inference  using  Markov  chain  Monte  Carlo  sampling  requires  a 
large  number  of  computations  of  the  source-receptor  relationship  to  be  undertaken.  In  con¬ 
sequence,  a  fast  and  efficient  technique  for  performing  computations  of  the  source-receptor 
relationship  (for  any  given  source  distribution  0)  is  required  for  the  rapid  sampling  of  the 
posterior  distribution.  To  this  purpose,  Keats  et  al  [14]  and  Yee  et  al  [17]  described  a 
computationally  efficient  methodology  for  determination  of  the  source-receptor  relationship 
using  an  adjoint  representation  for  this  relationship. 

Figure  2  illustrates  two  mathematically  equivalent  (dual)  representations  for  the  source- 
receptor  relationship  A4sr.  If  we  interpret  the  source  distribution  S  (which  we  encode 
as  0)  as  a  vector  in  V  (vector  space  of  source  functions),  then  we  can  construct  a  vector 
(function)  C*  which  is  dual  (or  conjugate)  to  S  belonging  to  the  dual  vector  space  V*  (or, 
conjugate  concentration  function  space).  This  vector  space  is  identified  to  be  the  space  of 
all  linear  functionals  C*  :  V  — >  M.  There  is  a  one-to-one  correspondence  between  the  vector 
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S  £  V  and  the  dual  vector  C*  £  V* ,  and  this  correspondence  can  be  defined  through  the 
scalar  product  ((C*\S))  that  pairs  S  with  C* . 3 

Similarly,  C  and  h  [see  Eq.  (9)]  can  be  interpreted  as  dual  (or  conjugate)  vectors  lying 
in  the  concentration  function  space  U  and  detector  function  space  U*,  respectively,  with 
U*  being  the  dual  vector  space  to  U  with  the  one-to-one  correspondence  between  ele¬ 
ments  of  these  two  spaces  defined  through  the  scalar  product  {{h\C))  with  C  £  U  and 
h  £  U* .  More  interestingly,  C  can  be  paired  with  its  dual  h  and  C*  can  be  paired  with 
its  dual  S  such  that  the  duality  relationship  ((C\h))  =  (( h\C ))  =  ((C*\S))  is  exactly  sat¬ 
isfied.  Now,  S  can  be  related  to  C  through  the  mapping  G  (viz.,  C  =  G(S)  =  GS) 
which  corresponds  to  the  direct  (source-oriented)  representation  of  the  source-receptor  re¬ 
lationship.  However,  a  mathematically  equivalent  representation  of  the  source-receptor 
relationship  can  be  formulated  by  relating  h  to  C*  through  the  adjoint  mapping  G*  (viz., 
C*  =  G*(h)  =  G*h )4  with  G*  explicitly  constructed  so  that  the  duality  relation  is  exactly 
satisfied:  ((C*\S)}  =  ((G*h\S))  =  ((h\GS))  =  {(h\C))  =  (( C\h )),  for  any  source  S  £  V  and  any 
receptor  h  £  U* .  This  is  the  dual  (receptor-oriented)  representation  for  the  source-receptor 
relationship. 

In  more  concrete  terms,  the  dual  (adjoint)  representation  of  the  source-receptor  relationship 
given  in  Eq.  (9)  can  be  expressed  explicitly  as  follows: 

C(xdi ,  tf ;  0)  =  [  dt  f  dx.  C(x,  t)h(x,  t\xd. ,  tf  )  =  (( C\h)){xdi ,  tf  ) 

3  Jto  Jv  3  3 

T 

=  J  dt'  J  dx'  C*  (x',  t'\xdi,t^.)S(x',  t') 

— oo  T> 

=  ((C*\S))(xdi,tf.),  (20) 

where  C'*(x',t/ |x^,^)  is  an  adjunct  (dual)  “concentration”  at  space-time  point  (x',f') 

(i) 

associated  with  the  sensor  concentration  data  at  location  xdi  and  time  tf] ' .  In  the  source- 

oriented  approach,  the  computation  of  Cj  using  Eqs  (9)  and  (1)  requires  the  determination 
of  C(x,t).  The  latter  quantity  requires  the  specification  of  the  kernel  function  p(x,  t|x',  t'). 
Unfortunately,  we  have  nothing  approaching  a  comprehensive  theory  of  turbulent  diffusion 
and,  as  a  consequence,  there  are  no  known  exact  solutions  for  p(x,  t\x',t')  for  the  case  of 
complex  turbulent  flows  (e.g.,  atmospheric  flows).  Approximations  for  the  kernel  function 
p(x,  f|x',f')  can  be  obtained  using  either  an  Eulerian  or  Lagrangian  description  of  atmo¬ 
spheric  diffusion. 

In  this  report,  we  will  use  a  first-order  Lagrangian  stochastic  (LS)  trajectory  simulation 
method  to  approximate  p(x,  f|x',f')  (and  thence,  C{x.,t)).  To  this  purpose,  we  consider 

3In  fact,  the  scalar  product  (((7*15)}  can  be  interpreted  as  an  explicit  mathematical  representation  for 
collection  of  all  linear  functionals  C*  :  V  — >  R  that  can  be  defined  on  V. 

4In  functional  analysis,  a  linear  operator  G*  :  U*  —>  V*  is  called  the  dual  or  pull-back  of  the  linear 
operator  G  :  V  -►  U  if  {(G*h\S))  =  ({h\GS)),  VS  £  V,  he  U*. 
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Thomson’s  [26]  model  for  inhomogeneous  Gaussian  turbulence,  where  the  forward  incre¬ 
ments  of  the  velocity  U(f)  =  ( Ui(t ))  and  position  X(t)  =  ( Xx(t ))  (i  =  1,2,3)  of  a  marked 
fluid  element  (particle)  are  given  by  the  following  stochastic  differential  equation:5 


dX(t)  =  U  (t)dt, 

dU{t)  =  a(X(*),U(t),t)<ft  +  (C0e(X(t),t))1/2  dW(t),  (21) 


where  e  is  the  turbulence  kinetic  energy  (TKE)  dissipation  rate,  Co  is  the  Kolmogorov 
“universal”  constant,  (ZW(i)  =  ( dWx(t ))  are  the  increments  of  a  vector- valued  (three- 
dimensional)  Wiener  process,  and  a  =  (a*)  is  the  drift  coefficient  vector  which  assumes 
the  following  form: 

=  -l{Coe)T-k}{uk  -  Uk)  +  (22) 

where 


Pe 


1  dTa  dUi  -  dUi 
+  -^L  +  Ul- 


2  dxx 

+ 


at 

V1 

2  b 


dxi 

ar„  ar. 


dt 


dxr 


+ 


(9.x, 


0xk 


(uj  Uj ) 


(23) 


Here,  C*  is  the  mean  Eulerian  velocity,  T^  =  ('iq  —  Ui)(uj  —  Uj)  is  the  Reynolds  stress 
tensor  (where  an  over  line  is  used  to  denote  an  ensemble  average) ,  and  Pe  is  the  background 
(Eulerian)  velocity  PDF  (which  is  implicitly  assumed  here  to  possess  a  Gaussian  form).  In 
this  source-oriented  approach,  marked  particles  with  initial  space-time  coordinates  (xj,tj) 
are  sampled  from  a  space-time  density  function  that  is  proportional  to  the  source  distribu¬ 
tion  5(x,  t),  and  the  forward-time  trajectories  of  these  tagged  fluid  particles  are  computed 
using  Eqs  (21)  to  (23)  for  t>  tt  (where  tt  is  the  initial  time  at  which  a  tagged  particle  was 
released  from  S(x,t)).  The  displacement  statistics  of  these  particles  are  used  to  determine 
C(x,f). 


For  source  reconstruction,  we  use  the  receptor-oriented  approach  described  by  Eq.  (20)  for 
the  efficient  calculation  of  Cj,  which  requires  the  determination  of  the  adjunct  (or  dual) 
concentration  C*  (x7,  t/|x(iJ,  tdj)  ■  To  this  purpose,  it  was  shown  by  Thomson  [26]  and  Flesch 
et  al  [27]  that  the  following  backward-time  Lagrangian  trajectory  simulation  model  is  the 
dual  to  the  forward-time  Lagrangian  trajectory  simulation  model  given  by  Eqs  (21)  to  (23): 

dXb{t')  =  U  b{t')dt', 

dUb{t')  =  a6(Xfe(t'),U b(t'),t')dt' +  (C0e(X.b(t'),t'))1/2  dW(t'),  (24) 

with 

*The  forward  increments  of  the  marked  particle  position  and  velocity  are  defined  as  dX(t)  =  ~K(t  +  dt)  — 
X(t)  and  d\J(t)  =  U(t  +  dt)  —  U (t),  respectively,  with  dt  >  0. 
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In  Eq.  (24),  dXb{t')  =  X.b{t')  -  X.b{t'  -  dt')  and  dUb{t')  =  U b{t')  -  U b{t'  -  dt ')  (, dt '  >  0)  are 
the  backward  increments  of  the  position  ~Kb{t')  and  velocity  U b(tr)  of  a  marked  particle  at 
time  t'  along  a  backward-time  trajectory.  In  this  receptor-oriented  (dual)  approach,  marked 
particles  with  final  space-time  coordinates  are  sampled  from  a  space-time  density  function 
that  is  proportional  to  the  sensor  response  function  /i(x;,  t'\xd,  td)  [at  the  detector  space- 
time  location  (xf/,  td)],  and  the  backward-time  trajectories  of  these  tagged  fluid  particles  are 
computed  using  Eqs  (24)  and  (25)  for  t'  <  tf  (where  tf  is  the  final  time  at  which  a  tagged 
particle  was  released  from  h(x.',  t'\x.d,  td)).  The  displacement  statistics  of  these  particles  are 
used  to  determine  C*(x',  t'\xd,  td)- 

Substituting  Eq.  (6)  into  Eq.  (9),  the  model  concentration  “seen”  by  the  detector  at  space- 
time  point  (x^ ,  tjj )  is  given  explicitly  as  follows: 

-  -  rmin(T,Tefc) 

Cj{Q)  =  C^di,tf-Q)  =  yjQk  C*{xs,k,ts\xdi,tf)dts.  (26) 

J  Jn 

It  should  be  noted  that  the  computation  of  Cj(Q)  can  be  obtained  for  any  source  distribu¬ 
tion  S  (encoded  as  the  parameter  vector  0),  without  having  to  re-compute  C*.  Indeed,  it 
should  be  emphasized  that  because  C*  does  not  depend  on  the  source  distribution,  it  can 
be  pre-calculated  using  the  backward-time  LS  trajectory  simulation  model  for  each  avail¬ 
able  detector  space-time  location,  and  this  pre-calculated  C*  can  be  used  in  Eq.  (26)  for  a 
computationally  efficient  determination  of  Cj(@)  (J  =  1,2, . . .  ,N)  required  for  the  rapid 
evaluation  of  the  likelihood  function  p(D|0,J). 

3.2  Markov  chain  Monte  Carlo  sampling 

All  the  information  arising  from  the  application  of  Bayesian  probability  theory  to  the 
problem  of  source  reconstruction  is  embodied  in  the  posterior  PDF  of  the  parameters 
e  =  (Ns,eNa)  eiw+1  [see  Eq.  (17)]  that  define  the  source  distribution  S.  The  posterior 
quantities  of  interest  are  expectation  values  with  a  generic  form  given  by  Eq.  (19),  which 
involves  an  integration  of  the  product  of  an  arbitrary  function  iP(Qns)  and  p(Ns,  9ns  |D,  I) 
(for  a  fixed  value  of  Ns)  over  a  subset  of  the  parameter  space  involving  0ns  6  M6Ars.  Un¬ 
fortunately,  in  potentially  high-dinrensional  spaces  ( Ns  S>  1),  a  numerical  integration  of 
Eq.  (19)  involving  the  evaluation  of  the  integrand  on  some  grid  in  the  8ns  -space  would  be 
prohibitively  expensive  (and  indeed,  even  for  the  case  Ns  =  1  implying  8ns  G  M6,  the  “curse 
of  dimensionality”  renders  the  quadrature  practically  infeasible). 

One  method  for  overcoming  this  “curse  of  dimensionality”  is  given  by  the  application  of 
Markov  chain  Monte  Carlo  (MCMC)  algorithms  for  posterior  sampling  (see  Gilks  et  al  [28] 
and  Gelrnan  et  al  [29]).  To  this  purpose,  Yee  [21]  described  the  formulation  of  a  reversible- 
jump  MCMC  (RJMCMC)  algorithm  applied  with  parallel  tempering  for  generating  samples 
from  the  posterior  distribution  p(0|D,7)  =  p(Ns,  0^s  |D,  I)  given  by  Eq.  (17).  As  a  con¬ 
sequence,  only  the  important  details  of  the  algorithm  that  are  required  to  understand  the 
following  Bayesian  analysis  of  field  trial  data  will  be  presented  here. 
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The  objective  of  MCMC  sampling  is  to  construct  an  auxiliary  Markov  chain  whose  sta¬ 
tionary  (or,  invariant)  distribution  is  the  posterior  distribution  p(0|D,7)  of  the  source 
parameters  0  =  ( Ns,0ns ).  The  difficulty  in  the  construction  of  the  Markov  chain  in  the 
current  application  resides  in  the  fact  that  number  of  sources  Ns  in  unknown  a  priori, 
so  the  dimension  of  the  hypothesis  (or,  parameter)  space  is  unknown  (viz.,  6ns  £  M6Ar,J 
with  Ns  unknown).  More  specifically,  with  reference  to  the  posterior  distribution  p(0|D,  7) 
given  by  Eq.  (17),  it  is  necessary  to  consider  (AiS)inax  —  AT,  +  1)  candidate  models  for  the 
source  distribution  S  and  associated  with  each  of  these  candidate  models  is  a  posterior  dis¬ 
tribution  p(Ns,  #ats|D,  7)  depending  upon  an  unknown  parameter  vector  0jvs  G  M.Ns  where 
Ns  G  {AT,  min,  AT,  m,n  +  1, . . . ,  AT,, max}  is  a  model  indicator  that  defines  the  dimension  ( 6NS ) 
of  the  hypothesis  (parameter)  space.  In  order  to  allow  changes  in  the  dimensionality  of  the 
model,  a  reversible-jump  MCMC  algorithm  is  used  to  construct  the  Markov  chain  for  0. 
The  formalization  of  RJMCMC  algorithms  for  dealing  with  variable  dimension  models  has 
been  described  by  Green  [30]. 

Consider  a  Markov  chain  with  “state”  vector  {©(*)}  =  {N^\  0  (t)}  (t  =  0, 1, 2, . . .)  that  is 
constructed  so  that  its  stationary  distribution  coincides  with  p(0|D,7)  given  by  Eq.  (17). 
After  convergence  of  the  Markov  chain  to  its  stationary  distribution,  the  samples  drawn  from 
this  chain  can  be  used  to  estimate  any  posterior  statistic  of  interest.  The  construction  of 
the  Markov  chain  uses  a  RJMCMC  algorithm  in  which  the  MCMC  moves  are  separated  into 
two  categories;  namely,  propagation  moves  which  do  not  change  the  dimensionality  of  the 
source  distribution  and  trans-dimensional  jump  moves  which  change  the  source  distribution 
by  ±1  discrete  source  (source  atom).  In  the  current  application,  the  trans-dimensional  jump 
move  changes  the  dimensionality  of  the  hypothesis  space  by  ±6  (as  each  source  atom  in 
Eq.  (6)  involves  six  degrees  of  freedom). 


For  the  dimension-conserving  moves,  we  partition  the  parameter  vector  9ns  (Ns  fixed)  as  fol¬ 
lows:  9Ns  =  (0\02)  where  01  =  (Qi,  Q2,  ■  ■  ■ ,Qns  )  €  M N°  and  92  =  (xS)i,  T^,  T*, . . .  ,xs,Nsr 
T^s ,Tgs)  G  M5JVs.  The  parameters  associated  with  91  are  linearly  related  to  the  model 
concentration  data,  as  is  evident  from  Eq.  (26).  For  these  parameters,  a  Gibbs  sampler  is 
used  for  the  update  (propagation)  move.  More  specifically,  the  Gibbs  sampler  updates  Qk  as 
a  direct  draw  from  the  univariate  full  conditional  posterior  distribution  p(Qk\Q-b,  92,  D,  7) 
where  9]_k  is  the  vector  01  with  its  fc-th  component  removed  (k  =  1,2 ,...,NS).  From 
Eq.  (17),  it  can  be  shown  that  the  full  conditional  posterior  distribution  for  Qk  assumes 
a  simple  Bernoulli-Gaussian  (truncated)  distribution  which  can  be  sampled  from  directly 
(see  Yee  [21]). 


The  remaining  parameters  in  92  (e.g.,  source  location,  source  on  time,  source  off  time) 
are  related  non-linearly  to  the  model  concentration.  In  consequence,  the  Gibbs  sampler 
cannot  be  used  to  update  these  parameters  owing  to  the  fact  that  the  full  conditional 
posterior  distribution  for  these  parameters  cannot  be  determined  analytically.  Furthermore, 
even  if  this  distribution  can  be  determined  analytically,  it  would  correspond  to  a  non¬ 
standard  probability  distribution  from  which  it  is  not  easy  (straightforward)  to  draw  a 
random  sample.  For  these  reasons,  the  Metropolis-Hastings  (M-H)  sampler  is  used  to  update 
the  parameters  in  92.  To  this  purpose,  we  update  92  G  02  (l  =  1, 2  . . .  ,  5 Ns)  by  randomly 
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sampling  a  new  candidate  Of  from  a  proposal  distribution  that  is  taken  to  be  a  mixture  of 
Gaussian  distributions,  each  with  mean  Of  and  different  variances  /3? .  The  variances  in  this 
mixture  of  Gaussian  distributions  are  chosen  typically  to  cover  several  orders  of  magnitude. 
For  the  M-H  steps  here,  we  used  a  mixture  consisting  of  seven  Gaussian  distributions  with 
the  base  standard  deviation  fli  for  a  particular  parameter  Of  taken  to  be  equal  to  0.01 
times  the  ‘length’  of  the  domain  of  definition  for  that  parameter  (assigned  using  the  prior 
distribution),  with  the  remaining  six  standard  deviations  chosen  to  be  equal  to  0.1,  0.2,  0.5, 
2,  5,  and  10  times  the  base  standard  deviation.  The  M-H  acceptance  probability  for  this 
update  move  for  Of  is  evaluated.  If  the  move  is  accepted,  the  new  candidate  Of  is  taken 
to  be  the  updated  value  for  this  parameter;  otherwise,  the  original  value  Of  becomes  the 
updated  value  for  this  parameter  (viz.,  the  Markov  chain  does  not  move  in  the  $p-subspace) . 

The  dimension-changing  moves  that  modify  the  source  distribution  by  ±1  source  atom 
(and  the  dimensionality  of  the  hypothesis  space  by  ±6)  are  provided  by:  (1)  a  creation 
move  C  that  results  in  the  addition  of  a  single  source  atom,  so  0'  =  (Ns  +  l,0/vs+i)  = 
C(0)  =  C((Ns,0ns))  where  0'  G  M1+6Afs+6;  and,  (2)  an  annihilation  move  C t  that  involves 
the  removal  of  a  single  existing  source  atom  from  the  current  source  distribution,  so  0r  = 
( Ns  —  1 , 6]\fs _ i )  =  Ct(0)  =  C^((Ns,6n J)  where  0'  G  M1+67Vs_6.  If  a  creation  move  C 
is  selected,  the  “coordinates”  of  the  new  source  atom  are  obtained  by  drawing  random 
samples  from  a  proposal  density  that  is  chosen  to  be  the  prior  density  for  each  coordinate. 
On  the  other  hand,  if  an  annihilation  (reverse)  move  C t  is  selected,  a  source  atom  in  the 
current  source  distribution  is  randomly  picked  and  removed.  The  acceptance  probability 
for  the  creation  move  C  or  the  annihilation  move  C 1  is  similar  to  that  for  a  M-H  algorithm 
involving  only  dimension-conserving  moves,  with  the  exception  of  a  term  involving  the 
Jacobian  of  the  transformation  C  or  C^.  From  this  viewpoint,  the  RJMCMC  is  simply  a 
generalized  M-H  sampler  for  the  hypothesis  space  that  includes  the  dimensionality  of  the 
space  as  an  additional  parameter  to  be  sampled.  Finally,  the  probabilities  for  creation 
and  annihilation  moves  are  chosen  as  follows:  pf?s  =  0  for  Ns  =  NSjmax  and  pfff  =  0  for 
Ns  =  Ns  min;  otherwise, 


Pcs 


vNa+1 
Pc  t 


1 

-  mm 
2 

1 

-  min 
2 


p(Ns  +  l\I)\ 

’  P(JV„|J)  /’ 

,  P(NS\I)  } 

'p(Ns  +  l\I)j- 


(27) 


To  summarize,  the  Markov  chain  consists  of  a  sequence  of  states  0^  (t  =  0,1,2,...) 
resulting  from  individual  updates  consisting  of  three  basic  moves:  (1)  involving  the 

creation  of  a  source  atom  at  a  random  location,  or  annihilation  of  an  existing  source  atom 
with  probabilities  pc  and  pc t,  respectively;  (2)  involving  updates  of  the  emission  rates 
of  the  source  atoms  using  Gibbs  sampling;  and,  (3)  M.%  involving  updates  of  the  location, 
source  on  time  and  source  off  time  of  the  source  atoms  using  M-H  sampling.  The  state 
vector  ©h^1)  of  the  Markov  chain  at  iteration  t  —  1  is  updated  to  the  state  vector  0^)  at 
iteration  t  using  the  following  procedure: 

1.  Specify  values  for  (IVs,min,  iVs,max,  Qmax,  Tmax,  t0,  J,p*)  which  define  p(Q\I). 
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2.  Choose  an  initial  state  for  the  Markov  chain  by  sampling  from  p(0|I). 

3.  For  i  6  { 1,2,...,  fUpper},  conduct  the  following  sequence  of  moves: 


©(*“!)  0*  ^  0**  ^  &(t),  (28) 

where  0*  and  0**  denote  some  intermediate  transition  states  between  iterations  t  —  1  and 
t. 

To  improve  the  “speed”  with  which  a  Markov  chain  traverses  the  hypothesis  space  (or, 
to  increase  the  mixing  rate  of  the  chain  in  the  hypothesis  space),  Yee  [21]  implemented  a 
form  of  parallel  tempering  based  on  a  Metropolis-coupled  MCMC  algorithm  described  by 
Geyer  [31].  In  this  approach,  r  Markov  chains  are  run  in  parallel,  each  with  a  different 
stationary  distribution.  These  chains  are  run  simultaneously,  but  occasionally  a  proposal 
is  made  to  swap  the  states  of  two  randomly  selected  chains.  In  consequence,  the  states  in 
the  “ladder”  of  Markov  chains  can  swap  positions  with  a  certain  acceptance  probability  as 
each  chain  equilibrates.  The  family  of  interrelated  stationary  distributions,  corresponding 
to  this  “ladder”  of  r  Markov  chains,  is  chosen  to  have  the  following  form: 

pi(0|D,/)ocp(0|7)pAi(D|0,/),  i  =  l,2,...,r,  (29) 

where  A,;  £  [0,1]  (i  =  1,2, ...  ,r)  is  an  increasing  sequence  (viz.,  A*  <  A j  for  i  <  j )  with 
Ai  =  0  and  Ar  =  1.  Note  that  each  distribution  in  Eq.  (29)  involves  a  tempering  parameter  A 
which  is  used  to  raise  the  likelihood  to  the  A  power  to  give  a  modified  posterior  proportional 
to  p(0|/)pA(D|0, 1).  When  A  =  0,  the  likelihood  function  is  switched  off  and  the  modified 
posterior  distribution  reduces  exactly  to  the  prior  distribution.  On  the  other  hand,  when 
A  =  1  the  modified  posterior  distribution  is  exactly  the  posterior  that  we  wish  to  sample 
from.  In  between  these  two  extremes,  with  A  £  (0,1),  the  effects  of  the  concentration 
data  D  are  introduced  gradually  through  the  modified  (or  “softened”)  likelihood  function 
pA(D|0,/). 

In  this  study,  rather  than  use  a  parallel  tempering  scheme,  we  employ  a  related  (and 
simpler)  simulated  annealing  scheme  to  facilitate  chain  mobility  in  the  hypothesis  space. 
In  this  scheme,  we  consider  an  ensemble  of  Nmem  (typically  between  50  and  200)  of  source 
distributions  (or,  source  molecules)  that  have  been  randomly  drawn  from  the  modified 
posterior  p^(0|D,7)  oc  p(0|/)pA(D|0, /).  These  samples  will  be  labelled  0fc(A),  with  A  £ 
[0,1]  ( k  =  1, 2, . . . ,  JVmem).  Note  that  po(@|D,7)  =  p(Q\I)  (prior  distribution  of  0)  and 
pi(@|D,/)  =  p(0|D,7)  (posterior  distribution  of  0).  In  this  framework,  it  is  useful  to 
interpret  the  tempering  parameter  A  as  an  inverse  temperature  parameter  T  (so,  A  =  1/T), 
with  A  £  [0,1]  implying  T  £  [1, oo] .  The  posterior  distribution  pi(0|D,7)  corresponds 
to  the  temperature  T  =  1,  whereas  the  modified  pa(0|D>7)  (A  £  [0,1))  corresponds  to 
“heating”  the  posterior  distribution  to  a  temperature  T  =  1/A  >  1  which  results  in  a 
flattening  of  the  distribution. 

When  the  stochastic  sampling  scheme  begins  and  A  =  0  (infinite  temperature) ,  we  randomly 
draw  Nmem  source  molecules  0^(0)  (k  =  1, 2  ... ,  Nmem)  from po(0|D,  I)  (prior  distribution); 


DRDC  Suffield  TR  2009-040 


15 


viz.,  ©fc(0)  ~  p(©|/).6  Given  an  ensemble  of  7Vmem  source  molecules  @fc(A)  that  has  achieved 
equilibrium  (at  temperature  T  =  1/A)  with  respect  to  the  modified  posterior  ^(€>113, /)/ 
an  ensemble  of  Nmem  source  molecules  @fc(A  +  5A)  that  is  consistent  with  Pa+<5a(©|E>,  I)  (at 
the  reduced  temperature  T  =  1  / (A  +  <5A),  <5A  >  0)  can  be  obtained  by  using  the  weighted 
resampling  method  (see  Gamerman  and  Lopes  [32])  applied  to  @^(A)  ( k  =  1,2, . . . ,  Nmem ) • 
To  this  purpose,  each  source  molecule  ©&( A)  is  assigned  an  importance  weight  Wk  as  follows: 


wk 


i>vw(8t(A)|D.J)  fc  1 

W(8t(A)|D,J)  \U  ’) 

(N mem  \  _1 

E  wA 


(30) 


for  k  =  1,2, .. . ,  Nmem .  Next,  a  resample  is  drawn  from  the  discrete  distribution  concen¬ 
trated  at  the  sample  of  source  molecules  {©/.'(A)}^!™  from  py(@|D,I)  with  respective 
weights  In  other  words,  the  resampling  step  here  involves  generating  a  new 

set  {©fc* }k*m=T  by  resampling  (with  replacement)  7Vmem  times  from  the  set  {@fc(A)}^Ym 
so  that  Pr{@fc*  =  0fc(A)}  =  w^.  This  resample  can  be  interpreted  as  an  ensemble  of 
source  molecules  &k*  (k*  =  1,2 ,lVmem)  that  is  drawn  from  the  modified  posterior 
Pa+<5a(©|D,  I). 

There  are  many  ways  to  implement  this  resampling  (with  replacement)  from  {©*.( A)}^jm, 
but  in  this  report  we  use  the  systematic  resampling  scheme  described  by  Kitagawa  [33]. 
This  scheme  requires  0(IVmem)  time  to  execute  and  minimizes  the  Monte  Carlo  vari¬ 
ation  in  the  resample.  The  basic  operations  for  this  scheme  can  be  described  briefly 
as  follows.  A  vector  (m,  ri2,  ■  ■  ■ ,  njvmem)  °f  copies  of  the  source  molecules  ©fc(A)  (k  = 
1,2...,  Nmem)  is  obtained  by  computing  a  vector  (pi,  p2, .  ■  ■ ,  PAfmem)  of  the  cumulative  sums 
of  Armem  (w\ ,  vj‘2 , . . .  •  w Nuiein ) ,  generating  a  random  draw  u  ~  ZL([0, 1]),  and  determining 

(k  =  1,2,...,  Nmem)  from 


k  —  \_Pk  T  u\  Ypk— i  T  •  k  —  2,3,...,  Nmem  1, 

m  =  [p!  +  u\,  njvmem  =  LpATmem-i  +u\,  (31) 

where  |_  •  J  denotes  the  “integer  part  of” .  After  resampling,  the  weights  for  each  member  of 
the  resample  are  reset  to  =  1/Nmem  (viz.,  an  equal  weight  is  assigned  to  each  member 
of  the  resample) .  After  each  resampling  step  to  give  an  ensemble  of  lVmem  source  molecules 
0fc(A  +  <5A)  that  is  in  equilibrium  at  temperature  T  =  1/(A  +  <5A)  (approximately  or  better) 
with  respect  to  Pa+sa(©|E>,  /),  we  apply  Nr  (typically  25)  Markov  chain  transitions  T\+§\  to 
each  member  of  the  ensemble  obtained  from  the  resampling  operation.  These  Markov  chain 
transitions  utilize  the  update  scheme  given  by  Eq.  (28)  and  leave  pa+5a(@|D>  I)  invariant. 

6The  prior  distribution  p(Q\I)  is  composed  of  standard  probability  distributions  for  which  independent 
sampling  is  easy. 

7This  simply  implies  that  {0fc(^)}fc/Tim  can  be  interpreted  as  an  ensemble  of  source  molecules  drawn 
from  Pa(6|D,  I). 
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An  annealing  schedule  for  A  G  [0, 1]  is  required  for  the  simulated  annealing.  In  this  report, 
we  applied  simulated  annealing  with  200  values  of  A  uniformly  spaced  in  the  interval  [0,  0.05] 
and  400  values  of  A  geometrically  spaced  in  the  interval  (0.05,1].  This  gentle  annealing 
schedule  allows  the  ensemble  of  Nmem  source  molecules  to  transition  slowly  through  a  series 
of  quasi-equilibrium  states  from  the  prior  distribution  (A  =  0,  or  infinite  temperature)  at  one 
end  of  the  annealing  schedule  to  the  posterior  distribution  (A  =  1,  or  unit  temperature)  at 
the  other  end  of  the  schedule.  When  A  =  1,  the  annealing  phase  is  complete  and  probabilistic 
exploration  of  the  hypothesis  space  proceeds  (for  each  of  the  lVmem  source  molecules  in  the 
ensemble)  in  accordance  to  the  scheme  summarized  in  Eq.  (28).  The  annealing  phase  of  the 
scheme,  corresponding  to  values  of  A  G  [0,1),  is  associated  with  the  burn-in  phase  of  the 
algorithm.  When  A  =  1,  the  MCMC  algorithm  has  reached  an  equilibrium,  at  which  point 
the  probabilistic  exploration  corresponding  to  the  sampling  from  the  posterior  distribution 
p(0|D,7)  begins.  These  samples  drawn  from  the  posterior  distribution  can  be  used  to 
make  inferences  about  all  characteristics  of  the  source  parameters  (e.g.,  posterior  means, 
variances,  HPD  intervals). 

Interestingly,  the  simulated  annealing  phase  of  the  proposed  scheme  can  be  used  to  esti¬ 
mate  the  normalization  constant  (or,  evidence)  Z  =  p(D|7)  for  the  posterior  distribution 
p(0|D,7)  [see  Eq.  (12)].  To  accomplish  this  objective,  we  use  an  approach  referred  to  as 
thermodynamic  integration  (see,  von  der  Linden  [34]),  which  has  its  origins  in  the  problem  of 
the  evaluation  of  partition  functions  in  statistical  mechanics  familiar  to  physicists.  Briefly, 
thermodynamic  integration  focusses  on  the  modified  evidence  Z( A)  defined  as  follows: 

Z(X)  =  J p(Q\I)p\B\Q,  I)  d&,  (32) 

which  can  be  seen  to  be  simply  the  normalization  constant  (evidence)  corresponding  to  the 
modified  posterior  pA(0|D,7).  Obviously,  Z( 0)  =  1  because  the  prior  p(Q\I )  is  normalized 
and  Z(l)  =  Z  is  the  desired  normalization  constant  (evidence)  for  the  posterior  p(0|D,7). 
Taking  the  logarithmic  derivative  of  Eq.  (32)  with  respect  to  A,  and  re-arranging  yields 

HZ(l)}  -  ln[Z(0)]  =  ln[Z(l)]  =  ['  <ln[p(D|0, 7)])A  dX,  (33) 

Jo 

where  (  •  }\  denotes  the  mathematical  expectation  operation  with  respect  to  the  modified 
posterior  pA(€>|Dj  7). 

Additionally,  the  computation  of  Z( A)  allows  the  determination  of  the  information  (gain) 
corresponding  to  the  receipt  of  the  concentration  data  D  and  the  updating  of  our  state  of 
knowledge  concerning  the  (unknown)  source  S  (encoded  as  0)  from  the  prior  distribution 
p(0|7)  to  the  posterior  distribution  p(0|D,7).  This  information  gain  (amount  of  useful 
information  about  0  embodied  in  D)  is  given  by  the  Kullback-Leibler  divergence  TAkl(A) 
at  A  =  1  [viz.,  by  TAkl(I)]  where  TAkl(A)  is  defined  as  follows  (Cover  and  Thomas  [35]): 

^kl(A)  =  /ln(^^y^)pA(©|D,7)d0 

=  X  j  ln[p(D|0,7)]pA(0|D,7)d0-ln[Z(A)]  j pA(0|D,7)70 
=  A(ln[p(D|0,7)])A-ln[Z(A)],  (34) 


DRDC  Suffield  TR  2009-040 


17 


on  using  the  fact  that  pa(@|D,/)  =  p(©|/)pA(D|0,  I)/Z(X).  The  Kullback-Leibler  diver¬ 
gence  defined  in  Eq.  (34)  for  A  =  1  is  simply  the  negative  of  the  entropy  (negentropy)  of 
the  posterior  relative  to  the  prior  and,  as  such,  is  the  information  gain  provided  by  the  re¬ 
ceipt  of  the  concentration  data  D.  More  specifically,  the  information  gain  “compresses”  the 
posterior  relative  to  the  prior  so  that  Dkl(1)  can  simply  be  interpreted  as  the  logarithm  of 
the  volumetric  factor  by  which  the  prior  has  been  compressed  to  become  the  posterior  (the 
greater  this  compression,  the  greater  is  the  information  gain  provided  by  the  concentration 
data). 

4  FUSION  Field  Trial  2007  (FFT-07) 


The  FUsing  Sensor  Information  from  Observing  Networks  (FUSION)  Field  Trial  2007 
(FFT-07)  was  conducted  to  obtain  research-grade  concentration  data  from  a  network  (or 
array)  of  fast-response  detectors  resulting  from  releases  of  a  passive  tracer  from  single  and 
multiple  sources.  The  scientific  objective  of  this  field  campaign  was  to  acquire  a  compre¬ 
hensive  meteorological  and  dispersion  dataset  that  can  be  used  to  validate  methodologies 
developed  for  source  reconstruction.  Details  of  the  instrumentation  deployed  and  the  exper¬ 
iments  conducted  in  FFT-07  are  given  in  Storwold  [36],  so  only  a  brief  summary  of  FFT-07 
will  be  presented  here.  In  particular,  only  the  relevant  details  of  the  experiments  that  are 
required  for  the  interpretation  of  the  results  in  this  report  are  emphasized. 

The  experiments  in  FFT-07  were  carried  out  in  September  2007  at  Tower  Grid  on  US  Army 
Dugway  Proving  Ground,  Utah  about  2  km  west  of  Camel  Back  Ridge  on  the  Great  Salt 
Lake  Desert.  The  terrain  was  flat,  uniform,  and  homogeneous  with  short  grass  interspersed 
with  low  shrubs  that  are  between  0.25  to  0.75  nr  in  height,  providing  an  upwind  fetch  that 
is  uniform  and  unobstructed  for  5  km  or  more  in  a  wide  sector.  The  test  elevation  site 
was  1330  nr  above  mean  sea  level,  and  the  terrain  gradually  rises  towards  the  southeast 
by  about  1  nr  over  horizontal  distances  of  about  2000  nr.  The  easterly  through  southerly 
drainage  flows  that  predominate  during  the  early  morning  hours  at  the  site  originate  on  the 
higher  terrain  to  the  southeast  and  are  channeled  by  Canrel  Back  Ridge. 

In  all  the  experiments,  the  tracer  gas  used  was  propylene  (CsHq)  which  was  chosen  be¬ 
cause  of  its  low  toxicity,  relatively  high  vapor  pressure  (938  kPa  at  21°C),  and  low  ion¬ 
ization  potential  (9.73  eV),  hence  giving  the  concentration  detectors  good  sensitivity.  The 
concentration  detectors  used  were  fast-response  digital  photo-ionization  (dPID)  detectors 
manufactured  commercially  by  Aurora  Scientific  Inc.  (Aurora,  Ontario,  Canada).  These 
detectors  give  a  frequency  response  of  50  Hz  with  a  sensitivity  of  about  0.025  parts  per  mil¬ 
lion  (ppnr)  by  volume  of  propylene.  The  detectors  were  calibrated  regularly  over  their  entire 
operating  range  using  a  specially  designed  calibration  unit  mounted  in  a  small  rack.  This 
unit  allowed  the  user  to  mix  gas  from  propylene  calibration  gas  cylinders  with  air  to  gener¬ 
ate  a  user-selectable  gas  concentration.  A  least-squares  regression  fit  of  the  calibration  data 
to  a  second-order  polynomial,  c  =  ao  +  aiU  +  02 U2,  provided  an  overall  calibration  for  each 
detector,  where  V  is  the  analog-to-digital  (A/D)  output  of  the  detector,  c  is  the  absolute 
propylene  concentration  (ppm),  and  a*  (i  =  0, 1,  2)  are  the  fitted  calibration  constants. 
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In  the  experiments,  a  plume  was  formed  in  the  atmospheric  surface  layer  by  releasing 
propylene  from  one  or  more  (up  to  a  maximum  of  four)  purpose-designed  gas  dissemination 
systems  (see  Chandler  [37]).  The  system  consisted  of  up  to  four  propylene  cylinders  con¬ 
nected  in  parallel  and  immersed  in  a  warm- water  bath.  A  regulator  was  connected  to  the 
outlet  of  the  cylinders  to  ensure  a  constant  downstream  pressure,  and  a  flexible  hose  was 
used  to  connect  the  regulator  to  the  inlet  fitting  of  the  mass  flow  controller.  This  controller 
was  used  to  set,  control,  and  measure  the  flow  through  the  dissemination  system.  The  con¬ 
troller  allowed  a  constant  gas  flow  rate  to  be  maintained  at  a  user-selected  reference  level 
between  10  and  250  1  min”1.  A  quick-release  connector  mated  the  outlet  of  the  mass  flow 
controller  to  the  dissemination  hose,  which  was  connected  to  the  base  of  the  disseminator, 
a  section  of  40  PVC  (polyvinyl  chloride)  pipe  1  nr  in  length  and  0.05  nr  in  diameter. 

The  network  (or  array)  of  concentration  detectors  used  in  FFT-07  is  shown  in  Figure  3.  A 
total  of  100  dPIDs  (indicated  by  the  filled  squares  in  Figure  3)  was  arranged  in  a  staggered 
configuration  consisting  of  10  rows  of  10  detectors.  The  rows  of  detectors  were  spaced  50  nr 
apart.  The  spacing  between  detectors  along  each  row  was  50  nr.  The  (x,  y)  local  coordinate 
system  that  will  be  used  here  is  shown  in  Figure  3.  The  direction  corresponding  to  the 
negative  x-axis  (referred  to  as  grid  north)  of  the  array  is  oriented  25°  west  of  north  to  take 
advantage  of  the  prevailing  wind  direction  at  the  test  site  in  the  sector  from  south-south- 
west  to  south-east.  Winds  incident  on  the  array  of  detectors  from  this  sector  resulted  in 
a  flat  and  homogeneous  upwind  fetch  of  more  than  10  km.  Note  that  x  is  the  lengthwise 
distance  from  the  windward  edge  (grid  south  edge)  of  the  array  and  y  is  the  spanwise 
distance  from  the  grid  west  edge  of  the  array.  The  overall  length  (along  the  x-direction) 
and  width  (along  the  y-direction)  of  the  detector  array  were  450  nr  and  475  nr,  respectively. 
Finally,  the  z  (vertical)  coordinate  is  defined  so  z  =  0  is  the  ground  surface. 

For  ease  of  reference,  the  rows  of  the  detector  array  will  be  numbered  from  1  to  10,  with 
row  1  forming  the  grid  north  edge  of  the  array  (at  x  =  —450  nr)  and  row  10  forming  the 
grid  south  edge  of  the  array  (at  x  =  0  nr).  Along  row  i  of  the  array  (z  =  1,2,...,  10),  the 
dPID  detectors  are  numbered  as  10(z  —  1)  +  j  (j  =  1,2,...,  10),  where  the  j  index  of  the 
detector  increases  from  grid  west  to  grid  east.  Consequently,  in  this  numbering  scheme  for 
detectors,  detector  number  91  [first  detector  ( j  =  1)  in  row  10  (i  =  10)]  is  located  at  the 
origin  of  the  local  (x,  y)  coordinate  system.  For  geo-referencing  purposes,  it  is  noted  that 
this  location  corresponds  to  the  following  geodetic  coordinates:  40.0923°  N  latitude  and 
112.9757°  W  longitude.  The  concentration  detectors  along  the  ten  sampling  lines  in  the 
array  were  placed  at  a  height,  Zd,  of  2.0  nr. 

In  the  experiments  used  for  the  current  analysis,  propylene  gas  was  released  continuously 
over  a  period  of  10  min  from  one  or  more  source  locations  (up  to  four)  at  a  height,  zs,  of 
2.0  nr.  The  four  source  locations  (labelled  1  to  4  in  Figure  3)  are  as  follows:  (1)  source 
1  is  at  (xs,ys)  =  (33.0,171.0)  nr;  source  2  is  at  (xs,ys)  =  (33.8,240.7)  nr;  source  3  is  at 
(xs,ys)  =  (30.0,312.9)  nr;  and,  source  4  is  at  (xs,ys)  =  (26.0,384.4)  nr.  Unfortunately,  for 
these  experiments,  only  the  mass  flow  controller  for  source  3  functioned  properly.  The  mass 
flow  controllers  for  sources  1,  2  and  4  failed  to  properly  regulate  the  flow  owing  to  the  fact 
that  the  impurities  in  the  low-grade  of  propylene  used  for  these  experiments  contaminated 
the  sensor  and  electronic  circuitry  in  these  controllers.  Consequently,  the  control  signals 


DRDC  Suffield  TR  2009-040 


19 


Source-Receptor  Configuration 


500  p 
450 
400 
350- 
300 
£  250- 
^200- 
150- 
100- 
50- 
0- 

--i 


•  4 

•  3 

•  2 

•  1 


-400  -300  -200  -100  0 

x  (m) 


100 


Figure  3:  A  schematic  diagram  of  the  geometry  of  the  network  (or,  array)  of  concentration 
detectors  (dPIDs)  used  in  FFT-07.  The  locations  of  the  100  detectors  (filled  squares)  and 
four  sources  (filled  circles)  are  shown. 


that  were  used  to  set  the  valve  to  regulate  the  flow  rate  in  these  controllers  failed  to  function 
correctly. 

Three-dimensional  (3-D)  sonic  anemometers  (R.  M.  Young,  Model  81000)  were  arranged 
on  three  32-m  lattice  towers  along  a  transect  parallel  to  the  x-axis  and  midline  of  the 
concentration  detector  array.  These  three  towers  were  located  at  the  center  of  the  detector 
array  at  (x,y)  =  (—225.0,237.5)  m  (or,  grid  center)  and  at  positions  750.0  m  upwind 
and  downwind  of  grid  center  at  (x,y)  =  (525.0,237.5)  m  and  (x,y)  =  (—975.0,237.5)  m, 
respectively.  Each  of  these  towers  was  instrumented  with  five  3-D  sonic  anemometers  at 
the  2-,  4-,  8-,  16-,  and  32-m  levels.  The  sonic  anemometer  data  were  recorded  at  10  Hz. 

An  estimate  of  the  momentum  roughness  length,  zq,  was  determined  from  mean  wind 
profiles  measured  under  near-neutral  stratification  (for  which  the  mean  wind  speed  variation 
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with  height  can  be  represented  by  a  semi-logarithmic  relation).  An  estimate  of  zo  =  1.3T0.2 
cm  was  obtained  from  the  3-D  sonic  anemometry  data  on  the  upwind  tower  at  (x,  y )  = 
(525.0,237.5)  m,  as  an  average  over  35  periods  (each  with  a  10-minute  sampling  time)  with 
the  magnitude  of  the  Obukhov  length  L  from  the  3-D  sonic  anemometer  at  2-m  height 
exceeding  500  nr  (viz.,  \L\  >  500  nr  at  z  =  2  nr). 

Table  1:  Summary  of  turbulence  statistics  measured  with  a  3-D  sonic  anemometer  at  the 
2-m  level  on  the  upwind  tower.  Here,  S2  is  the  horizontal  mean  wind  speed  at  the  2-m 
level,  a  is  the  mean  wind  direction,  u*  is  the  friction  velocity,  L  is  the  Obukhov  length, 
and  <Ji  (i  =  u,v,w)  are  the  standard  deviations  of  the  fluctuating  wind  velocity  in  the  three 
coordinate  directions  (alongwind,  crosswind,  and  vertical,  respectively). 


Trial 

S2 

(nr  s-1) 

a 

(deg) 

u* 

(nr  s-1) 

L 

(nr) 

®uf 

(-) 

(Tv  / 

(-) 

^  e 

1 

e 

* 

I 

2.60 

142.3 

0.418 

10.7 

2.18 

1.62 

1.17 

II 

3.00 

138.9 

0.183 

15.4 

2.36 

1.72 

1.23 

III 

3.61 

155.1 

0.282 

-27.3 

2.33 

1.86 

1.10 

The  mean  wind  and  turbulence  statistics  data  for  three  experiments  used  to  test  the  source 
reconstruction  methodology  are  given  in  Table  1.  Here,  the  wind  components  are  calculated 
in  a  rotated  coordinate  system  with  the  alongwind  (u),  crosswind  (v),  and  vertical  (w)  com¬ 
ponents,  aligned  along  the  x'-,  y1-,  and  z-directions,  respectively.  The  standard  deviations 
of  the  fluctuating  wind  in  these  three  coordinate  directions  are  denoted  ou,  crv,  and  aw , 

1  /  c\ 

respectively.  Measurements  of  the  mean  horizontal  wind  speed  S2  =  ( u 2  +u2)  and  mean 
wind  direction  a,  obtained  from  the  3-D  sonic  anemometer  at  the  2-m  level  on  the  upwind 
tower,  are  summarized  in  Table  1  for  these  three  experiments.  The  mean  wind  direction 
corresponding  to  a  wind  from  the  + x-axis  direction  (see  Figure  3)  will  be  denoted  in  the 
usual  compass  convention  as  a  =  155°  (and  is  associated  with  a  wind  direction  from  grid 
south).  This  mean  wind  direction  corresponds  to  normal  incidence  on  the  detector  array 
(viz.,  the  direction  is  perpendicular  to  the  sampling  lines  of  detectors  along  the  y-axis) . 
For  simplicity,  a'  will  denote  the  deviation  of  the  mean  wind  direction  from  a  =  155°  (or, 
normal  incidence).  The  Obukhov  length  L  was  estimated  as 

r  ufT 

L  = - - - , 

Kgw'T' 


where  T  =  T  +  T'  is  the  sonic  temperature  (which  can  be  decomposed  into  a  mean  value 
T  and  fluctuation  therefrom),  u *  is  the  friction  velocity  at  the  surface,  kv  ~  0.4  is  the 
von  Karrnan  constant,  and  g  is  the  acceleration  due  to  gravity.  The  friction  velocity  was 
estimated  as  follows: 


u*  =  \{u'w')2  +  {v'w')2  1^4 . 


(36) 
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5  Application  to  FFT-07  data 


In  this  section,  we  apply  the  source  reconstruction  algorithm  to  the  three  dispersion  data 
sets  identified  in  Table  1.  Trials  I,  II  and  III  correspond  to  one-  three-  and  four-source 
examples. 

In  all  three  examples,  the  backward-time  LS  model  given  by  Eqs  (24)  and  (25)  was  used 
to  determine  Cj(@)  from  Eq.  (26).  This  model  was  applied  to  short-range  dispersion  in 
the  atmospheric  surface  layer  over  a  level  and  unobstructed  terrain.  The  mean  wind  and 
turbulence  statistics  for  the  three  examples  were  assumed  to  be  horizontally  homogeneous 
and  stationary,  with  the  relevant  flow  statistics  for  these  examples  summarized  in  Table  1. 
The  backward-time  LS  model  was  applied  with  the  Kolmogorov  constant  Co  =  4.8,  a 
value  which  was  recommended  by  Wilson  et  al  [38]  from  a  calibration  of  the  model  against 
concentration  data  obtained  from  Project  Prairie  Grass. 

The  wind  flow  and  turbulence  statistics  parameterization  for  the  LS  model  are  prescribed 
based  on  standard  surface-layer  relationships  from  Monin-Obukhov  theory  and  summarized 
in  the  Appendix  of  Yee  [21].  These  relationships  have  been  modified  slightly  so  as  to  be 
consistent  with  the  measured  values  of  au/u *,  a v/u*  and  aw/u *  (at  the  surface)  compiled 
in  Table  1  for  the  three  trials.  Finally,  the  three  surface-layer  parameters  zo,  u *  and  L 
required  for  the  wind  statistics  parameterization  have  been  either  reported  in  main  text  of 
Section  4  or  summarized  in  Table  1. 

All  the  examples  in  this  report  involve  continuously  emitting  sources  and,  as  a  consequence, 
— >  — oo  and  Tefc  — >  oo  in  Eqs  (6)  and  (26)  [. k  =  1,2,...  ,NS].  As  a  consequence,  the 
relevant  source  parameters  are  as  follows:  0  =  (Ns,xsi,Qi, . . . ,  xstvs ,  Q Na ) •  Furthermore, 
it  is  assumed  that  the  height  of  the  sources  above  ground  level  (zs  =  2.0  nr)  is  known  a 
priori ,  so  the  only  unknown  location  parameters  are  (xs,ys)  of  the  sources  in  the  horizontal 
plane.  In  view  of  these  assumptions,  the  adjunct  concentration  C*(xs|xd)  in  Eq.  (26)  can 
be  pre-calculated  for  one  detector  position  (for  the  known  source  height  zs)  in  any  given 
trial,  with  the  adjunct  concentration  C*  (considered  as  a  function  of  xs)  at  all  other  detector 
positions  obtained  by  a  simple  translation  of  C'*(xs|x(i)  in  the  horizontal  (x,y)  plane. 

5.1  Trial  I:  one-source  example 

For  this  example,  the  source  corresponding  to  source  3  (see  Figure  3)  was  turned  on.  The 
(constant)  emission  rate  from  this  source  was  Q  =  2.8  g  s'1.  The  source-detector  config¬ 
uration  for  this  example  is  depicted  in  Figure  4.  Only  six  detectors  (shown  by  the  filled 
blue  squares  in  Figure  4)  in  the  array  (viz.,  detector  numbers  25,  36,  46,  56,  66,  and  77) 
were  used  for  the  source  reconstruction.  For  this  example,  the  mean  wind  direction  was 
a'  =  155°  —  a  =  12.7°  away  from  normal  incidence  to  the  detector  array. 

It  is  assumed  that  the  number  of  sources  is  known  a  priori  (viz.,  this  knowledge  formed 
part  of  the  background  information  I)  for  this  example.  To  this  end,  Ns  min  =  -YSjinax  =  l.8 

8Note  that  the  value  of  the  hyperparameter  p*  that  characterizes  the  binomial  prior  for  Ns  in  this  example 
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Figure  4:  In  Trial  I,  one  source  located  upwind  of  the  array  of  detectors  was  turned  on.  The 
solid  dot  shows  the  location  of  this  source.  Squares  show  the  location  of  the  detectors  in 
the  array:  open  and  filled  squares  indicate  that  the  detector  at  the  given  location  is  missing 
and  present,  respectively,  in  the  array.  A  filled  blue  square  marks  the  detectors  that  were 
used  for  the  source  reconstruction. 

The  other  parameters  needed  to  define  the  prior  distribution  p(Q\I)  for  this  example  are 
chosen  as  follows:  Qmax  =  100  g  s_1;  7  =  1.0;  and  V  =  [0, 100]  x  [0, 500]  providing  the  prior 
bounds  on  the  source  location  in  the  (x,  y)-plane.  Recall  that  the  (x,  y )  coordinate  system 
used  here  is  chosen  as  shown  in  Figure  4. 

After  the  initial  simulated  annealing  phase  was  completed  using  Nmem  =  50  members  of 
an  ensemble  of  source  distributions  (originally,  randomly  drawn  from  the  prior  distribution 
p(Q\I )  for  A  =  0),  the  MCMC  algorithm  was  run  for  an  additional  1,000  iterations  for  each 
member  of  the  ensemble,  giving  a  total  of  50,000  samples  of  source  distributions  (encoded  as 
0)  obtained  from  the  posterior  distribution  p(@|D,I)  during  the  probabilistic  exploration 
phase  (A  =  1)  of  the  algorithm.  The  marginal  distributions  for  the  source  location  ( xs,ys ) 
and  the  emission  rate  Q  =  qs  obtained  from  these  50,000  samples  are  shown  in  Figure  5. 
The  actual  source  location  was  (. xs ,  ys)  =  (30.0,  312.9)  m  and  the  emission  rate  was  Q  =  2.8 
g  s_1.  The  analysis  of  the  samples  gave  the  following  estimates  for  the  source  parameters 
expressed  as  a  posterior  mean  value  and  standard  deviation  (s.d.)  and  the  boundary  (lower 

is  irrelevant,  owing  to  the  fact  that  p(Ns\I)  =  1  ( Ns  G  [IVs.min, iYs,max])  for  As, min  =  Af,,max  =  1. 
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Figure  5:  The  marginal  posterior  histograms  for  the  source  location  (xs,ys)  a nd  for  the 
emission  rate  qs  obtained  for  Trial  I.  The  vertical  solid  lines  mark  the  true  values  of  the 
source  parameters  and  the  vertical  dashed  lines  mark  the  best  estimates  of  the  source 
parameters  obtained  from  the  posterior  mean. 


and  upper)  of  the  95%  HPD  interval:  xs  =  19.9  ±  11.6  m  (0.0,40.5)  nr;  ys  =  311.2  ±  1.3  nr 
(309.0,313.6)  nr;  and,  qs  =  3.9  ±  0.6  g  s~x  (2.7, 5.1)  g  s_1.  Although  the  information 
embodied  in  the  six  concentration  detectors  allowed  the  crosswind  position  ys  of  the  source 
to  be  well  determined,  this  information  only  weakly  constrained  the  alongwind  position  xs 
of  the  source  which  is  estimated  with  a  large  uncertainty.  As  a  consequence,  there  is  also  a 
large  uncertainty  in  the  determination  of  the  source  emission  rate.  Indeed,  the  information 
gain  (provided  by  the  concentration  data)  for  this  example  was  determined  to  be  only 
D kl(1)  =  10.2  natural  units  (nits)  [see  Eq.  (34)]. 

5.2  Trial  II:  three-source  example 

For  Trial  II,  three  sources  were  turned  on  upwind  of  the  detector  array  as  shown  in  Figure  6. 
These  sources  correspond  to  sources  2,  3  and  4  shown  in  Figure  3.  The  emission  rate  Q  from 
source  3  was  known  to  be  3.8  g  s-1,  whereas  the  emission  rates  for  sources  2  and  4  were  not 
known  in  this  example  owing  to  the  fact  that  the  mass  flow  controllers  for  these  two  sources 
malfunctioned  because  of  their  contamination  by  the  impurities  in  the  propylene  cylinders 
connected  to  them.  The  source  reconstruction  for  this  example  will  be  attempted  using  37 
detectors  in  the  array,  which  have  been  marked  using  a  filled  blue  square  in  Figure  6.  In 
this  trial,  the  mean  wind  direction  had  an  obliquity  angle  with  respect  to  the  normal  to  the 
detector  array  of  a'  =  155°  —  a  =  16.1°. 
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Figure  6:  In  Trial  II,  three  sources  located  upwind  of  the  array  of  detectors  were  turned 
on.  The  solid  dots  show  the  locations  of  the  three  sources.  Squares  show  the  location  of 
the  detectors  in  the  array:  open  and  filled  squares  indicate  that  the  detector  at  the  given 
location  is  missing  and  present,  respectively,  in  the  array.  A  filled  blue  square  marks  the 
detectors  that  were  used  for  the  source  reconstruction. 

For  this  example,  we  assume  that  the  number  of  sources  Ns  is  unknown  a  priori.  Further¬ 
more,  it  is  assumed  that  an  exact  upper  bound  on  the  maximum  number  of  possible  sources 
is  available.  In  consequence,  for  our  specification  of  the  binomial  prior  for  Ns,  we  choose 
IVSimin  =  1  and  Ns, max  =  4,  with  p*  =  1/3. 9  As  in  the  first  example,  the  hyperparameters 
Q max  and  7  that  define  the  Bernoulli-uniform  prior  for  the  emission  rate  Q,  are  initialized 
to  100.0  g  s'1  and  0.25,  respectively.  The  prior  bounds  on  the  (xs,ys)  location  of  each  of 
the  unknown  sources  are  given  by  V  =  [0, 100]  x  [0,500]. 

We  applied  our  source  reconstruction  algorithm  to  this  example.  The  simulated  annealing 
phase  of  the  algorithm  was  initiated  (at  A  =  0)  with  A^mem  =  50  different  source  distributions 
S  (encoded  as  0)  drawn  randomly  from  the  prior  distribution  p(Q\I).  A  sequence  of 
modified  posterior  distributions  with  A  E  (0, 1]  and  the  associated  annealing  schedule  as 
described  in  Section  3.2  was  used  in  the  simulated  annealing  phase.  After  the  termination 
of  this  phase  with  A  =  1,  a  further  1000  iterations  of  the  RJMCMC  procedure  were  applied 

9This  choice  for  p(N3\I)  implies  that  the  expected  number  of  sources  in  the  domain  is  (NB)  =  A8>min  + 
(As, max  —  AS!min )p*  =  2.  In  consequence,  the  prior  distribution  for  As  favors  the  wrong  choice  for  the 
number  of  sources  in  this  example. 
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Figure  7:  The  posterior  distribution  for  the  number  of  sources,  p(Ns)  =  p(Ns\D,I),  for 
Trial  II  estimated  using  50,000  samples  obtained  from  the  probabilistic  exploration  phase 
of  the  stochastic  sampling  algorithm. 

to  each  of  the  -/Vmem  =  50  members  of  the  ensemble  of  source  distributions  to  give  50,000 
samples  of  source  distributions  drawn  from  the  posterior  distribution  p(@|D,  /)  during  the 
probabilistic  exploration  phase  of  the  algorithm. 

Figure  7  displays  the  posterior  probability  distribution  for  the  number  of  sources  p(Ns)  = 
p(-/Vs|D,I)  for  this  example.  The  most  probable  number  of  sources  is  Ns  =  3,  and  this 
coincides  with  the  correct  number  of  sources.  In  consequence,  our  maximum  a  posteriori 
estimate  for  the  number  of  discrete  sources  present  in  Trial  II  is  Ns  =  3  [cf.  Eq.  (18)],  and 
this  choice  is  favoured  with  a  probability  of  about  0.82.  It  is  important  to  emphasize  that 
samples  with  Ns  =  4  discrete  sources  were  also  drawn  from  the  posterior  distribution  during 
the  probabilistic  exploration  phase  of  the  algorithm.  However,  it  turns  out  that  with  the 
prior  specification  7  =  PrjQfc  >  0}  =  0.25  (k  =  NStm[n,  NS:U1[n  +  1, . . . ,  lVSimax),  many  of  the 
samples  obtained  with  Ns  =  4  had  one  of  the  discrete  sources  turned  off  (so,  the  emission 
rate  Q  =  0  for  this  source).  For  the  case  of  a  four-source  sample,  but  with  one  of  the  sources 
turned  off  (with  the  result  that  this  source  does  not  contribute  to  the  model  concentration 
“seen”  by  the  detectors  in  the  array),  we  classify  the  source  distribution  sample  here  as 
having  three  discrete  sources  rather  than  four.  This  convention  for  the  determination  of 
the  number  of  discrete  sources,  associated  with  a  sample  of  a  source  distribution  model, 
was  used  in  the  construction  of  Figure  7. 

Given  the  fact  that  our  best  estimate  for  the  number  of  sources  is  Ns  =  3,  we  can  now 
estimate  the  source  parameters  ( xs,ys,qs )  corresponding  to  each  of  these  three  discrete 
sources.  Towards  this  objective,  we  first  extract  all  samples  of  source  distributions  drawn 
from  the  posterior  distribution  p(@|D,I)  having  exactly  three  sources  (viz.,  we  choose 
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Figure  8:  Inference  of  the  discrete  source  parameters  obtained  from  samples  drawn  from 
the  posterior  distribution  p(@|D,  I )  having  exactly  three  discrete  sources,  (a)  Samples  from 
the  posterior  distribution  having  three  discrete  sources  projected  onto  the  (xs,ys)  subspace. 
(b,c,d)  Histograms  for  the  three  parameters,  namely  alongwind  location  xs,  crosswind  lo¬ 
cation  ys,  and  emission  rate  qs  that  characterize  sources  2,  3,  and  4  (cf.  Figure  3).  In 
each  frame,  the  solid  vertical  line  indicates  the  true  value  of  the  parameter  (if  known )  and 
the  dashed  vertical  line  corresponds  to  the  best  estimate  of  the  parameter  obtained  as  the 
posterior  mean  of  the  marginal  posterior  distribution  for  the  parameter. 


all  sample  indices  t  such  that  has  a  first  component  that  satisfies  Ns  =  3  for  t  = 
1,2, .. .  ,50000).  Figure  8(a)  displays  samples  corresponding  to  Ns  =  3,  projected  onto  the 
(xs,ys)  subspace.  From  this  figure,  we  can  identify  three  clusters  of  points  that  determine 
the  locations  of  the  three  discrete  sources.  Observe  that  the  concentration  data  D  constrain 
the  crosswind  positions  ys  of  the  three  sources,  but  there  is  considerable  uncertainty  in  the 
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determination  of  the  alongwind  positions  xs  of  the  sources  (which  is  especially  acute  for  the 
source  located  at  the  largest  value  of  ys). 


Table  2:  The  posterior  mean,  posterior  standard  deviation,  and  lower  and  upper  bounds 
of  the  95%  HPD  interval  of  the  parameters  xs ^  (m),  Vs,k  (m),  and  qs ^  (g  s~1)  for  k  = 
1,  2,  and  3  calculated  from  samples  of  source  distribution  models  with  Ns  =  3  (the  latter 
corresponding  to  the  most  probable  number  of  sources  in  the  domain  as  inferred  from 
Figure  7). 


Parameter 

Mean 

Standard  Deviation 

95%  HPD 

Actual 

k  = 

1 

xs  (m) 

30.7 

9.3 

(11.0,47.5) 

33.8 

Vs  (m) 

239.0 

2.7 

(233.5,243.5) 

240.7 

Qs  (g  s'1) 

7.4 

0.7 

(5.96,8.81) 

- 

k  = 

2 

xs  (nr) 

30.7 

6.2 

(18.6,41.6) 

30.0 

Vs  (nr) 

313.1 

2.3 

(308.5,317.2) 

312.9 

Qs  (g  s_1) 

4.1 

0.5 

(3.1, 5.1) 

3.8 

k  = 

3 

xs  (nr) 

37.5 

19.8 

(0.8,73.0) 

26.0 

Vs  (nr) 

388.3 

5.7 

(379.0,399.3) 

384.4 

Qs  (g  s_1) 

6.5 

1.0 

(4.8, 8.9) 

- 

An  important  issue  in  the  interpretation  of  the  results  is  the  identifiability  problem  that 
arises  owing  to  the  fact  that  the  posterior  distribution  of  0  is  invariant  under  a  reordering 
(or,  relabelling)  of  the  identifiers  used  for  each  discrete  source  in  the  source  distribution 
model  [cf.  Eq.  (17)].  More  specifically,  it  can  be  seen  that  p(0|D,J)  is  invariant  under  a 
permutation  of  the  discrete  source  identifier  k.  However,  from  Figure  8(a),  it  is  evident  that 
we  can  uniquely  identify  the  discrete  sources  in  a  source  distribution  model  if  we  impose 
an  ordering  constraint  on  the  ys ^-locations  of  the  three  sources  so  that  yS)\  <  ys ^  <  ys, 3- 
Furthermore,  comparing  Figures  3  and  8(a),  it  can  be  seen  that  after  relabelling  the  discrete 
sources  of  each  source  distribution  model  in  accordance  to  this  ordering  constraint,  the 
sources  at  ySt\,  yS) 2  and  ys, 3  are  associated  with  sources  2,  3  and  4,  respectively.  Figures  8(b), 
(c)  and  (d)  display  the  histograms  of  the  alongwind  position  xs,  crosswind  position  ys  and 
emission  rate  qs  of  the  three  identified  sources,  using  the  same  label  for  the  sources  as 
provided  in  Figure  3.  The  posterior  mean  and  standard  deviation,  as  well  as  the  95% 
HPD  interval  of  the  discrete  source  parameters  for  each  of  the  three  identified  sources  are 
summarized  in  Table  2.  Note  that  sources  2  (k  =  1)  and  4  (k  =  3),  whose  emission  rates 
were  not  regulated  owing  to  the  failure  of  the  mass  flow  controllers  connected  to  these 
sources,  appear  to  have  higher  emission  rates  than  source  3  ( k  =  2)  whose  emission  rate 
was  properly  regulated.  Indeed,  the  measured  emission  rate  from  source  3  for  Trial  II  was 
Q  =  3.8  g  s-1,  which  compares  well  with  the  posterior  mean  estimate  for  this  emission  rate 
given  in  Table  2  (for  k  =  2).  Finally,  the  information  gain  provided  by  the  concentration 
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Figure  9:  In  Trial  III,  four  sources  located  upwind  of  the  array  of  detectors  were  turned 
on.  The  solid  dots  show  the  locations  of  the  four  sources.  Squares  show  the  location  of 
the  detectors  in  the  array:  open  and  filled  squares  indicate  that  the  detector  at  the  given 
location  is  missing  and  present,  respectively,  in  the  array.  A  filled  blue  square  marks  the 
detectors  that  were  used  for  the  source  reconstruction  (Case  1). 


data  for  this  example  was  found  to  be  Dkl{  1)  =  35.2  nits. 

5.3  Trial  III:  four-source  example 

5.3.1  Case  1 :  62  detectors 

This  example  involves  four  continuously  emitting  sources.  In  case  1,  the  detectors  used  for 
the  source  reconstruction  algorithm  are  shown  in  Figure  9.  As  can  be  seen  from  this  figure, 
this  case  involves  the  use  of  62  detectors  in  the  array  for  source  inversion.  All  the  detec¬ 
tors  in  the  array  that  measured  a  significantly  non-zero  mean  concentration  were  used  for 
the  reconstruction,  as  well  as  a  number  of  detectors  for  which  the  measured  mean  concen¬ 
tration  was  nominally  zero  (viz.,  over  the  sampling  time  for  the  trial,  the  instantaneous 
concentration  did  not  exceed  the  detection  threshold  concentration).  In  this  example, 
the  mean  wind  direction  was  normally  incident  to  the  detector  array;  more  specifically, 
a'  =  155°  —  a  =  —0.1°  away  from  normal  incidence  to  the  detector  array  (cf.  Table  1). 


DRDC  Suffield  TR  2009-040 


29 


Instantaneous  estimates  of  N 

S 


Probability  distribution  of  N ' 


Figure  10:  Trace  plot  (top)  of  the  number  of  discrete  sources  Ns  in  the  source  distribution 
model  samples  drawn  from  p(0|D,J)  during  the  probabilistic  exploration  phase  of  the 
stochastic  sampling  algorithm  for  Trial  III  (62  detectors  used  for  source  reconstruction), 
and  the  corresponding  posterior  distribution  for  the  number  of  sources,  p(Ns )  =  p(Ars|D,  I). 
In  the  trace  plot,  the  number  of  samples  displayed  has  been  decimated  by  a  factor  of  10. 


The  proposed  stochastic  sampling  algorithm  randomly  initializes  all  unknown  source  pa¬ 
rameters  in  accordance  to  the  prior  distribution  p(@|D).  In  this  example  (unlike  the  case 
dealt  with  in  Trial  II  above) ,  it  is  assumed  that  we  do  not  have  a  prior  knowledge  of  the  max¬ 
imum  number  of  available  sources  used  in  FFT-07.  As  a  consequence,  we  choose  IVS  m jn  =  1 
and  JV<  maY  =  8,  with  p*  =  1/7  in  the  specification  of  p(Ns\I)  implying  that  the  expected 
number  of  sources  is  (Ns)  =  2.  In  consequence,  our  initial  specification  for  the  prior  dis¬ 
tribution  of  Ns  favors  the  wrong  choice  for  the  actual  number  of  sources.  The  remaining 
hyperparameters  defining  p(@|D)  are  chosen  as  follows:  7  =  0.25  and  Qmax  =  100.0  g  s_1 
for  specification  of  the  prior  for  the  emission  rate  Q ;  and,  V  =  [0, 100]  x  [0, 500]  nr  which 
is  used  to  define  the  prior  bounds  for  the  location  ( xs,ys )  of  any  source.  An  ensemble  of 
Amem  =  50  members  of  source  distribution  models  0  were  drawn  from  p(@|D)  and  used 
for  the  simulated  annealing  phase  of  the  MCMC  algorithm.  After  A  =  1  was  achieved, 
1000  further  iterations  of  the  RJMCMC  algorithm  were  applied  to  each  of  these  source 
distribution  model  members  during  the  probabilistic  exploration  phase  of  the  algorithm  to 
give  50,000  samples  of  source  distribution  models  drawn  from  the  posterior  distribution 
P(B|D,I). 
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Figure  10  (top)  shows  a  trace  plot  for  the  number  of  discrete  sources  in  a  source  distribution 
model  against  the  sample  (or,  iteration)  number.  From  this  plot,  it  can  be  seen  that  after 
the  initial  burn-in  phase  of  the  stochastic  sampling  algorithm  (simulated  annealing  phase), 
the  samples  of  source  distribution  models  drawn  from  p(@|D,J)  during  the  probabilistic 
exploration  phase  generally  mixes  well  over  Ns.  In  particular,  annihilation  moves  for  models 
from  Ns  =  4  to  3  do  not  occur.  However,  dimension-changing  moves  involving  transitions 
from  Ns  =  4  to  5  (and,  vice-versa),  as  well  as  higher-order  transitions  (e.g.,  from  Ns  =  6 
to  7  and  its  reverse)  are  seen  to  occur.  Even  so,  it  is  seen  that  the  transitions  to  large 
values  of  Ns  (e.g.,  Ns  ~  lYs,max)  are  rare  and  short-lived.  Figure  10  (bottom)  displays  the 
marginal  posterior  distribution  p(Ns )  =  p(A^s|D,/)  for  the  number  of  sources.  Note  that 
the  most  probable  number  of  sources  for  Trial  III  is  4  ( Ns  =  4),  which  is  favoured  with  a 
probability  of  about  0.6.  The  most  probable  value  for  Ns  in  this  case  coincides  with  the 
correct  number  of  sources  Ns  =  4.  This  result  is  obtained  in  spite  of  the  fact  that  the 
prior  distribution  for  Ns  was  initialized  with  p*  =  1/7  (with  an  a  priori  expected  number 
of  sources  ( Ns )  =  2),  implying  that  the  algorithm  is  not  sensitive  to  this  hyperparameter. 
The  information  embodied  in  the  concentration  data  was  sufficient  to  move  the  stochastic 
simulations  towards  the  more  complex  model  with  Ns  =  4  source  atoms. 

Figure  11(a)  displays  samples  of  all  source  distribution  models  drawn  from  the  posterior 
distribution  p(0|D,/)  with  Ns  =  4.  Note  that  there  are  four  clusters  of  points,  with  the 
centroids  of  each  of  these  clusters  coinciding  (approximately  or  better)  with  the  true  location 
of  the  four  sources  (see  Figure  3).  It  is  evident  that  the  concentration  data  D  constrain 
the  ys-locations  of  the  four  sources,  but  the  xs-locations  of  these  sources  are  subject  to 
greater  uncertainty.  An  examination  of  Figure  11(a)  suggests  that  identifiability  of  the 
sources  can  be  obtained  by  ordering  the  discrete  sources  of  each  source  distribution  model 
on  the  ys-coordinate.  More  specifically,  the  label  k  for  each  discrete  source  in  a  source 
distribution  model  sample  are  reordered  so  that  the  ys ^.-locations  of  the  discrete  sources 
verify  ySj i  <  ys $  <  ys, 3  <  UsA-  With  this  relabelling  of  discrete  sources,  the  index  k  for 
a  source  follows  the  labelling  of  sources  given  by  Figure  3.  Figure  11(b)  shows  traces  of 
the  source  parameters  associated  with  each  discrete  source  of  the  source  distribution  model 
samples  (with  Ns  =  4),  after  the  relabelling  of  the  discrete  sources.  As  can  be  seen  from 
this  figure,  the  ys-locations  of  the  discrete  sources  are  well  separated,  but  the  xs-locations 
and  emission  rates  qs  of  the  discrete  sources  are  very  similar  and  significantly  overlap  one 
another. 

Figure  12  shows  the  marginal  posterior  distribution  (histogram)  of  the  parameters  xs,  ys 
and  qs  for  each  of  the  four  discrete  sources  identified  in  Figure  11.  The  posterior  mean  and 
standard  deviation,  as  well  as  the  lower  and  upper  bounds  for  the  95%  HPD  interval,  for 
each  of  these  four  identified  discrete  sources  are  summarized  in  Table  3.  For  this  example,  it 
is  seen  that  generally  the  estimates  for  the  source  parameters  are  quite  good  and,  certainly, 
the  true  values  of  the  parameters  (when  these  are  known)  lie  within  the  stated  errors. 

An  examination  of  Table  3  shows  that  the  emission  rates  for  the  three  unregulated  sources 
(viz.,  sources  1,2,  and  4)  generally  are  larger  than  that  for  the  regulated  source  (viz.,  source 
3).  More  specifically,  measurements  of  the  emission  rate  from  source  3  for  Trial  III  yielded 
Q  =  3.8  g  s-1,  which  agrees  with  the  estimated  value  given  in  Table  3  for  this  source  of 


DRDC  Suffield  TR  2009-040 


31 


Samples  Drawn  from  Posterior 


40  60 

x  (m) 


100 


(a)  Sample  density  plot 


60 


Iteration  N  umber  x  i  o4 

(b)  Trace  plots  of  source  parameters 

Figure  11:  ( a )  Density  plot  consisting  of  samples  of  source  distribution  models  obtained  for 
Ns  =  4  projected  onto  the  (xs,ys)  subspace,  (b)  Trace  plots  of  source  parameter  estimates 
against  sample  (or,  iteration)  number,  after  relabelling  of  discrete  sources  (+  (blue):  k  =  1; 
o  (black):  k  =  2;  □  (green):  k  =  3;  A  (red):  k  =  4).  In  the  trace  plots,  the  number  of 
samples  has  been  decimated  by  a  factor  of  100. 

qs  =  4.1  ±  0.4  g  s-1.  As  in  the  previous  examples,  the  localization  of  the  source  in  the 
x^-direction  (alongwind)  is  generally  poorer  than  that  in  the  ys-direction  (crosswind).  It 
can  be  seen  that  the  crosswind  locations  of  the  sources  can  be  determined  with  an  accuracy 
of  about  ±1  m  (standard  deviation),  but  the  alongwind  locations  of  the  sources  can  only 


32 


DRDC  Suffield  TR  2009-040 


Source  1 


Source  2 


v  (m) 


q.  fg  O 


x  (m) 


V  (m) 


(a) 


(b) 


Source  3  Source  4 


x  (m)  3;  (m)  (m)  3;  (m) 


?s(gS  ') 


3000 

2500 

2000 

1500 

1000 


?s(gs  ') 


(c) 


(d) 


Figure  12:  Inference  of  the  discrete  source  parameters  obtained  from  samples  drawn  from  the 
posterior  distribution  p(@|D,I)  having  exactly  four  discrete  sources.  (a,b,c,d)  Histograms 
for  the  three  parameters,  namely  alongwind  location  xs,  crosswind  location  ys,  and  emission 
rate  qs  that  characterize  sources  1,  2,  3,  and  4  (cf.  Figure  3).  In  each  frame,  the  solid 
vertical  line  indicates  the  true  value  of  the  parameter  (if  known)  and  the  dashed  vertical 
line  corresponds  to  the  best  estimate  of  the  parameter  obtained  as  the  posterior  mean  of 
the  marginal  posterior  distribution  for  the  parameter. 


be  determined  with  an  accuracy  of  about  ±5  nr  (standard  deviation).  This  appears  to  be 
true  in  this  example,  except  for  source  2  whose  alongwind  location  is  determined  to  an 
accuracy  of  about  ±2  m  (standard  deviation) .  The  reason  for  the  increased  accuracy  in  the 
determination  of  the  alongwind  location  of  source  2  stems  from  the  fact  that  only  detector 
96  (along  line  10)  -  see  Figure  9  -  measured  a  significantly  non-zero  mean  concentration  and 
this  concentration  was  due  solely  to  the  emission  from  source  2.  The  information  content 
embodied  in  this  detector  allowed  the  alongwind  location  of  source  2  to  be  inferred  with 
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Table  3:  The  posterior  mean,  posterior  standard  deviation,  and  lower  and  upper  bounds  of 
the  95%  HPD  interval  of  the  parameters  xs ^  (in),  ys,k  (m),  and  qs ^  (g  s1)  for  k  =  1,  2, 
3,  and  4  calculated  from  samples  of  source  distribution  models  with  Ns  =  4  for  Case  1  (the 
latter  corresponding  to  the  most  probable  number  of  sources  in  the  domain  as  inferred  from 
Figure  10). 


Parameter 

Mean 

Standard  Deviation 

95%  HPD 

Actual 

k 

=  1 

xs  (nr) 

34.0 

6.0 

(22.8,45.9) 

33.0 

Vs  (nr) 

170.9 

0.8 

(169.3,172.5) 

171.0 

Qs  (g  s_1) 

8.2 

0.6 

(7.0, 9.4) 

- 

k 

=  2 

xs  (nr) 

33.8 

1.6 

(30.5,37.0) 

33.8 

Vs  (nr) 

240.5 

0.4 

(239.8,241.5) 

240.7 

Qs  (g  s'1) 

7.1 

0.6 

(5.9, 8.4) 

- 

k 

=  3 

xs  (nr) 

23.7 

5.2 

(13.1,32.1) 

30.0 

Vs  (nr) 

313.4 

0.9 

(311.8,315.2) 

312.9 

Qs  (g  s_1) 

4.1 

0.4 

(3.3, 4.9) 

3.8 

k 

=  4 

xs  (nr) 

25.6 

4.3 

(17.5,33.8) 

26.0 

Vs  (nr) 

384.5 

0.7 

(383.2,385.8) 

384.4 

Qs  (g  s_1) 

6.1 

0.5 

(5.2,  7.0) 

- 

greater  accuracy  than  the  other  three  sources.  Finally,  for  this  example,  the  information  gain 
obtained  from  the  concentration  data  D  was  found  to  be  Djcl(  1)  =  52.5  nits,  implying  that 
the  information  contained  in  the  concentration  data  allowed  the  “posterior  volume”  of  the 
hypothesis  space  (volume  of  hypothesis  space  of  reasonably  large  plausibility  after  receipt 
of  the  concentration  data)  to  decrease  by  a  factor  of  exp(DA'L(l))  ~  6.3  x  1022  relative  to 
the  “prior  volume”  of  the  hypothesis  space  (volume  of  hypothesis  space  of  reasonably  large 
plausibility  before  the  receipt  of  the  concentration  data). 

5.3.2  Case  2:  45  detectors 

In  case  2  of  the  example  of  four  continuously  emitting  sources  (Trial  III),  we  consider  the  use 
of  45  detectors  (shown  in  Figure  13)  in  the  array  for  source  reconstruction.  In  contrast  to 
Case  1,  all  the  detectors  on  sampling  lines  8,  9  and  10  (viz.,  the  three  sampling  lines  closest 
to  the  emitting  sources)  have  been  removed.  The  stochastic  sampling  algorithm  was  applied 
to  the  concentration  data  obtained  from  the  remaining  45  detectors,  using  exactly  the  same 
hyperparameters  for  the  prior  distribution  p(Q\I)  as  described  previously  for  Case  1.  The 
simulated  annealing  phase  was  applied  to  an  ensemble  of  Nmem  =  50  members  of  source 
distribution  models  0  drawn  fromp(0|J).  After  the  termination  of  the  simulated  annealing 
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Figure  13:  In  Trial  III,  four  sources  located  upwind  of  the  array  of  detectors  were  turned 
on.  The  solid  dots  show  the  locations  of  the  four  sources.  Squares  show  the  location  of 
the  detectors  in  the  array:  open  and  filled  squares  indicate  that  the  detector  at  the  given 
location  is  missing  and  present,  respectively,  in  the  array.  A  filled  blue  square  marks  the 
detectors  that  were  used  for  the  source  reconstruction  (Case  2). 


phase,  the  RJMCMC  algorithm  was  run  for  1000  iterations  on  the  Nmem  =  50  members  of 
the  ensemble  (associated  with  the  inverse  temperature  A  =  1)  to  generate  50,000  samples 
of  source  distribution  models  drawn  from  the  posterior  distribution  p(@|D,/)  (during  the 
probabilistic  exploration  phase). 

Figure  14  exhibits  a  histogram  of  the  number  of  sources  Ns  for  this  case,  obtained  from 
the  50,000  samples  of  source  distribution  models  drawn  from  p(@|D,/).  The  RJMCMC 
simulations  settle  in  a  distribution  which  slightly  favors  Ns  =  5  sources.  However,  the 
hypotheses  Ns  =  4  and  6  are  equally  (approximately  or  better)  probable,  albeit  with  a 
smaller  probability  than  the  hypothesis  Ns  =  5.  The  results  of  Figure  14  would  suggest 
that  the  most  probable  number  of  source  in  Trial  III  is  Ns  =  5  (which  turns  out  to  be 
incorrect). 

Next,  let  us  extract  all  samples  of  source  distribution  models  having  exactly  Ns  =  5  discrete 
sources.  Figure  15(a)  displays  these  samples  projected  onto  the  (xs,ys)  plane.  There  are 
four  prominent  “streaks”  roughly  parallel  to  the  xs-axis,  and  the  narrowness  of  these  streaks 
usefully  constrain  the  locations  of  four  discrete  sources  in  the  ys-direction.  The  length  of 
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Figure  14:  The  posterior  distribution  for  the  number  of  sources,  p(Ns )  =  p(iVs|D,7),  for 
Trial  III  (Case  2)  estimated  using  50,000  samples  obtained  from  the  probabilistic  exploration 
phase  of  the  stochastic  sampling  algorithm. 


the  streaks  imply  that  the  xs-locations  of  these  sources  are  only  weakly  constrained  by  the 
concentration  data  D  (viz.,  there  is  significant  uncertainty  in  the  x^-locations  of  the  four 
identified  sources).  In  addition  to  the  four  streaks  of  points,  there  is  a  diffuse  cloud  of  points 
in  the  (xs,  ys)-plane  associated  with  the  location  of  the  fifth  discrete  source  in  the  source 
distribution  model  samples.  Note  that  these  points  appear  to  be  randomly  scattered  in  the 
(xs,  ys)-plane,  implying  that  the  concentration  data  D  do  not  contain  information  that  can 
usefully  constrain  the  location  of  this  fifth  discrete  source. 

Figure  15(b)  exhibits  the  marginal  posterior  distributions  for  the  source  parameters  corre¬ 
sponding  to  the  samples  of  source  distribution  models  with  Ns  =  5  discrete  sources.  Note 
that  the  histogram  for  ys  exhibits  four  modes  whose  locations  coincide  (approximately  or 
better)  with  the  ys -locations  of  the  four  sources  present  in  Trial  III.  The  histogram  for  xs 
has  a  single  mode  at  xs  ~  25  m  which  corresponds  roughly  to  the  xs-locations  of  the  four 
sources,  but  this  marginal  posterior  distribution  is  very  broad  indicating  a  large  uncertainty 
in  the  determination  of  the  alongwind  locations  of  the  sources.  Finally,  the  histogram  for 
qs  is  bimodal.  There  are  broad  modes  at  qs  ~  0  g  s^1  and  at  ~  7  g  s”1.  The  latter  mode  is 
associated  with  the  four  emitting  sources  identified  in  the  histogram  for  ys.  The  mode  with 
emission  rate  near  zero  is  associated  with  the  fifth  discrete  source,  and  it  is  seen  that  the 
emission  rate  from  this  source  is  small.  Indeed,  the  concentration  contributed  by  this  fifth 
source  is  generally  smaller  than  the  uncertainty  in  the  specification  of  the  concentration 
data  D.  This  simply  implies  that  this  source  is  not  usefully  constrained  by  the  available 
concentration  data  for  this  case,  and  may  simply  correspond  to  an  extraneous  source. 

Even  though  Figure  14  suggests  that  the  most  probable  number  of  sources  is  Ns  =  5,  the 
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Figure  15:  (a)  Density  plot  consisting  of  samples  of  source  distribution  models  obtained  for 
Ns  =  5  projected  onto  the  (xs,ys)  subspace,  (b)  Histogram  of  the  alongwind  location  xs, 
crosswind  location  ys,  and  emission  rate  qs  of  the  discrete  sources  obtained  from  samples  of 
source  distribution  models  with  Ns  =  5.  The  vertical  line(s)  in  each  panel  marks  the  true 
value  of  the  source  parameter  (if  known). 


analysis  of  the  samples  of  source  distribution  models  with  Ns  =  5  discrete  sources  (see 
Figure  15)  suggests  that  only  four  of  these  sources  are  usefully  constrained  by  the  con¬ 
centration  data  D.  Of  course,  there  may  have  been  one  or  more  additional  sources  than 
those  identified  in  Figure  15,  but  if  the  concentration  data  do  not  constrain  these  additional 
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Figure  16:  Inference  of  the  discrete  source  parameters  obtained  from  samples  drawn  from  the 
posterior  distribution  p(0|D,/)  having  exactly  four  discrete  sources.  (a,b,c,d)  Histograms 
for  the  three  parameters,  namely  alongwind  location  xs,  crosswind  location  ys,  and  emission 
rate  qs  that  characterize  sources  1,  2,  3,  and  4  (cf.  Figure  3).  In  each  frame,  the  solid 
vertical  line  indicates  the  true  value  of  the  parameter  (if  known)  and  the  dashed  vertical 
line  corresponds  to  the  best  estimate  of  the  parameter  obtained  as  the  posterior  mean  of 
the  marginal  posterior  distribution  for  the  parameter. 


sources,  then  no  useful  inferences  can  be  made  about  them.  As  a  consequence,  it  is  useful  to 
extract  all  samples  of  source  distribution  models  with  Ns  =  4  discrete  sources,  and  use  this 
information  to  infer  the  source  parameters  for  the  four  sources  identified  in  Figure  15.  To¬ 
wards  this  objective,  Figure  16  displays  the  histograms  of  the  source  parameters  {xs,  ys,  gs} 
constructed  from  all  samples  with  Ns  =  4  [after  relabelling  the  discrete  sources  for  each 
sample  so  that  their  ysy  locations  are  ordered  as  ysy  <  ysg  <  ys, 3  <  ys, 4  —  an  ordering 
constraint  that  is  obvious  to  impose  after  a  perusal  of  Figure  15(b)]. 
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Table  4:  The  posterior  mean,  posterior  standard  deviation,  and  lower  and  upper  bounds  of 
the  95%  HPD  interval  of  the  parameters  xs ^  (m),  ys^  (m),  and  qs ^  (g  s ~l)  for  k  =  1,  2,  3, 
and  4  calculated  from  samples  of  source  distribution  models  with  Ns  =  4  for  Case  2. 


Parameter 

Mean 

Standard  Deviation 

95%  HPD 

Actual 

k 

=  1 

xs  (m) 

46.7 

23.7 

(3.0,90.2) 

33.0 

Vs  (m) 

170.3 

2.6 

(165.2,175.3) 

171.0 

Qs  (g  s_1) 

8.5 

1.0 

(6.6,10.6) 

- 

k 

=  2 

xs  (m) 

35.3 

22.9 

(0.0,80.7) 

33.8 

Vs  (m) 

239.4 

2.9 

(233.6,244.8) 

240.7 

Qs  (g  s  1 ) 

6.8 

0.9 

(4.0, 8.6) 

- 

k 

=  3 

x s  (nr) 

50.1 

26.4 

(6.4,98.8) 

30.0 

Vs  (m) 

312.2 

3.8 

(304.5,319.2) 

312.9 

Qs  (g  s"1) 

4.5 

0.8 

(3.1, 6.1) 

3.8 

k 

=  4 

xs  (m) 

20.6 

10.3 

(0.4,38.8) 

26.0 

Vs  (m) 

385.6 

2.2 

(381.1,389.8) 

384.4 

Qs  (g  s_1) 

5.8 

0.6 

(4.6,  7.0) 

- 

Table  4  summarizes  the  posterior  mean  and  standard  deviation,  as  well  as  the  lower  and 
upper  bounds  for  the  95%  HPD  interval  of  the  source  parameters,  for  each  of  the  four 
discrete  sources.  The  estimated  values  for  the  crosswind  locations  of  the  four  sources  agree 
well  with  the  true  crosswind  locations,  albeit  the  uncertainties  in  these  inferred  locations 
are  now  about  three  times  as  large  as  those  for  Case  1.  The  estimates  of  the  alongwind 
locations  of  three  of  the  four  sources  are  quite  good,  although  the  uncertainties  in  these 
estimates  are  large.  An  examination  of  Figure  16(c)  and  Table  4  for  k  =  3  indicates  that  the 
available  concentration  data  D  for  Case  2  does  not  constrain  the  a^-location  for  Source  3. 
Indeed,  from  Figure  16(c),  it  can  be  seen  that  the  marginal  posterior  distribution  of  xs  for 
Source  3  is  almost  a  uniform  distribution  over  its  domain  of  definition  (and,  hence,  almost 
identical  to  the  marginal  prior  distribution  of  xs  for  Source  3).  The  posterior  mean  value 
of  qs  for  Source  3  overestimates  the  true  value  of  the  emission  rate  ( Q  =  3.8  g  s_1).  Even 
so,  the  true  emission  rate  here  is  contained  within  the  one  standard  deviation  uncertainty 
interval  for  this  parameter.  Finally,  for  Case  2,  the  information  gain  provided  by  the 
concentration  data  D  was  found  to  be  Dkl{  1)  =  39.2  nits  for  Case  2  (which  is  13  nits  less 
information  gain  than  provided  by  the  concentration  data  for  Case  1).  Viewed  in  another 
way,  the  concentration  data  available  on  the  last  three  sampling  lines  (which  were  used  for 
the  source  reconstruction  for  Case  1,  but  not  for  Case  2)  provided  13  nits  of  additional 
information,  and  this  information  allowed  the  posterior  “volume”  in  the  hypothesis  space 
for  Case  1  to  be  reduced  by  a  factor  of  exp(13)  ~  4.4  x  105  relative  to  the  posterior  “volume” 
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in  the  hypothesis  space  for  Case  2. 


6  Conclusions 


In  the  report,  we  have  developed  and  tested  an  innovative  Bayesian  method  for  source 
reconstruction  for  the  difficult  case  when  the  number  of  sources  is  unknown  a  priori.  More 
specifically,  Bayesian  probability  theory  was  applied  to  formulate  a  posterior  distribution  for 
the  number  of  sources  and  for  the  parameters  that  characterize  each  of  these  sources  (e.g., 
location,  emission  rate,  time  of  release,  etc.).  The  evaluation  of  this  posterior  distribution  in 
order  to  extract  its  features  of  interest  is  realized  using  an  efficient  computational  technology 
described  in  Yee  [21].  The  computational  algorithm  uses  simulated  annealing  for  the  burn-in 
phase  and  a  reversible-jump  MCMC  method  for  the  probabilistic  exploration  phase. 

The  source  reconstruction  methodology  has  been  successfully  validated  against  a  real  dis¬ 
persion  experiment  involving  various  combinations  of  multiple  source  releases  (FFT-07), 
with  measurements  of  the  resulting  concentration  field  obtained  from  an  array  of  100  detec¬ 
tors.  In  particular,  three  different  trials  were  used  to  test  the  methodology:  namely,  Trial  I 
was  used  to  test  the  source  reconstruction  methodology  for  the  case  in  which  the  number  of 
sources  is  known  a  priori ;  Trial  II  was  used  to  test  the  source  reconstruction  methodology 
for  the  case  in  which  the  number  of  sources  is  unknown  a  priori ,  but  the  maximum  possible 
number  of  sources  is  known  (four  in  this  case);  and,  Trial  III  was  used  to  test  the  source 
reconstruction  methodology  for  the  case  in  which  the  number  of  sources  is  unknown  a  pri¬ 
ori,  and  the  maximum  number  of  possible  sources  is  not  available.  For  Trial  III,  the  effect 
of  reducing  the  number  of  detectors  on  the  quality  of  the  inference  was  investigated.  These 
various  examples  illustrate  the  effectiveness  of  the  proposed  source  reconstruction  method¬ 
ology  and  demonstrate  the  reliable  determination  of  the  number  of  sources  and  estimation 
of  the  source  parameters  (along  with  the  associated  uncertainties)  corresponding  to  each  of 
the  identified  sources. 

There  are  several  possible  extensions  to  this  work.  A  number  of  the  field  experiments  con¬ 
ducted  in  FFT-07  involved  the  presence  of  a  spurious  additive  background  in  the  concentra¬ 
tion  measured  on  the  detector  array.  It  would  be  useful  to  extend  the  source  reconstruction 
methodology  to  treat  an  unrecognized  spurious  source  of  signal.  This  may  require  the  de¬ 
velopment  of  a  rational  procedure  within  the  Bayesian  framework  that  would  allow  the 
proper  separation  of  the  true  concentration  signal  from  an  underlying  background  (whose 
source  may  be  unknown).  The  various  contributions  to  the  noise  term  in  our  model  of  the 
mean  concentration  observations  were  simply  lumped  together  and  it  was  assumed  that 
a  noise  scale  parameter  (or  standard  deviation)  associated  with  the  aggregate  error  was 
known.  Of  course,  apparent  inconsistencies  in  the  inference  may  arise  from  an  incorrect 
estimate  of  the  uncertainty  of  this  aggregate  error.  It  would  be  useful  to  generalize  the 
methodology  to  treat  the  possible  uncertainty  in  the  specification  of  the  standard  deviation 
for  the  aggregate  error. 
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